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I .  INTRODUCTION 

Ever  since  the  generalized  network  formulation  for  aperture  problems 
was  given  in  terms  of  the  method  of  moments  [1],  solutions  for  particular 
problems  have  been  obtained  using  particular  subsections.  For  example, 
rectangular  patches  were  used  for  rectangular  apertures  [2],  and  annular 
subsections  were  used  for  annular  apertures  [3] .  Babinet's  principle  plus 
the  wire-grid  model  of  the  complementary  conducting  plate  have  been  used 
for  arbitrarily-shaped  apertures  [4].  This  approach  is  often  satisfactory 
for  far-field  quantities  and  transmission  coefficients,  but  is  not  appropri¬ 
ate  for  computing  near-field  quantities.  This  is  because  there  are  diffi¬ 
culties  in  relating  computed  wire  currents  to  equivalent  surface  magnetic 
currents.  Also,  the  accuracy  of  the  wire-grid  approximation  can  be  ques¬ 
tioned  on  theoretical  grounds. 

In  this  report,  the  problem  of  electromagnetic  transmission  through 
an  arbitrarily-shaped  aperture  in  an  infinite  conducting  screen  of  zero 
thickness  is  investigated  using  triangular  patches  to  model  the  aperture. 

The  method  of  solution  is,  in  general,  a  specialization  of  that  for  bodies 
of  arbitrary  shape  by  Rao  [5].  In  the  formulation,  the  equivalence  principle 
and  image  theory  [6]  are  used  to  derive  an  integral  equation  for  the  equi¬ 
valent  magnetic  currents.  The  moment  method  [7,8]  is  used  to  metricize 
this  integral  equation.  The  expansion  functions  are  chosen  to  be  local 
position  vectors  inside  each  triangular  patch. 

Extensions  of  the  basic  problem  are  also  given.  One  extension 
is  two  half  spaces  with  different  media.  Another  is  a  lossy  dielectric 
window  covering  the  aperture.  Computer  programs  are  written  and  numerical 
results  for  the  magnetic  currents,  transmission  cross  section  patterns  and 
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transmission  coefficients  are  given  for  several  sample  cases. 


II.  STATEMENT  OF  THE  PROTOTYPE  PROBLEM 

The  problem  configuration  to  be  considered  is  shown  in  Fig.  1. 

An  infinite  conducting  screen  with  an  arbitrarily-shaped  aperture  covers 
the  entire  xy-plane.  The  excitation  of  this  aperture  is  an  arbitrarily- 
polarized  plane  wave  incident  from  the  region  z  >  0  at  an  angle  0^  to  the 
z-axis.  The  quantities  to  be  computed  are  the  equivalent  magnetic  current 
distribution  and  the  transmission  characteristics  of  the  aperture. 

As  described  in  [1],  we  use  the  equivalence  principle  and  image 
theory  to  obtain  equivalent  situations  for  both  regions.  The  solution  is 
expressed  in  terms  of  the  equivalent  magnetic  current  M  =  E  x  z  in  the 
aperture.  To  compute  M,  we  use  a  linear  expansion  of  basis  functions  ^ 
and  moment  methods  to  evaluate  the  coefficients.  Hence,  we  have  to 
determine  a  generalized  admittance  matrix  and  an  excitation  vector.  To 
predict  the  transmission  characteristics,  we  need  a  measurement  vector. 
Since,  the  incident  field  is  a  plane  were,  the  excitation  vector  is  of 
the  same  form  as  the  measurement  vector. 


III.  FUNDAMENTAL  FORMULATION 

Refer  to  the  generalized  network  formulation  for  aperture  problems 

[1].  Define  M  =  J  V  M  over  the  aperture  region.  Then 
n 


(1) 


where 


[Y*l  -  (y‘'i  -  [y"*)  .  t  <- 


2[y'»,  .  21<  -  W,. 


(2) 


3 


Y 


Fig,  1.  Prototype  problem  configuration. 


V 


Hence, 


In  free  space,  the  magnetic  field  produced  by  a  source  M  is 


H(M  )  =  -  JwF  -V$ 


where  F  and  $  are  the  electric  vector  potential  and  the  magnetic  scalar 
n  n 


potential  related  to  M  as  follows  [6] 


4  ■  M„  •  0(k,t,r')ds' 


n  Airy 


m  G(k,r,r*)ds' 

Tl  —  — 


m  =  V  •  M 
n  ju)  — n 


where  the  free  space  Green's  function  is 


G(k,r,r')  = 


,-Jk|r-r’: 


Hence,  the  element  Y  in  the  admittance  matrix  is 
ran 


y®  +  y'’  =  4  <-  W  ,  H^®(M  )> 
mn  ran  m  t  — n 


4  ^  ,  (Ju^  +  V$^)ds’ 

4  |j  <J«s„  •  4  -  V  •  V'**' 


Ajo)  (W  •  F  +  w  $  )ds’ 
-m  -m  m  n 
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IV.  TRIANGULAR  PATCHES  AND  BASIS  FUNCTIONS 

Different  approaches  to  model  apertures  with  some  simple,  or  highly 
symmetrical  shapes  (e.g.  slot,  rectangle,  circle)  have  been  developed. 

Here  we  use  triangular  patches  for  the  sake  of  being  able  to  conform  closely 
to  arbitrarily-shaped  apertures.  There  are  other  advantages:  First,  this 
triangular  patch  scheme  is  easily  inputed  to  the  computer,  since  the  vertices 
can  be  independently  specified.  Second,  it  also  provides  the  flexibility  of 
having  greater  patch  densities  on  those  portions  of  the  aperture  where  more 
resolution  is  desired,  e.g.  when  we  are  concerned  about  the  edge  effect. 

The  presence  of  derivatives  on  the  magnetic  current  and  on  the 
scalar  magnetic  potential  suggests  that  we  have  to  be  careful  in  selecting 
the  expansion  functions  and  testing  procedures  in  the  method  of  moments. 

As  Rao  did  for  a  scattering  body  [5], we  choose  a  set  of  basis  functions 
which  yield  a  continuous  magnetic  current  and  a  piecewise  constant  mag¬ 
netic  charge  representation. 

Assume  that  a  suitable  trlangulatlon  defined  by  an  appropriate 

set  of  patches,  edges,  vertices,  and  boundary  edges,  such  as  shown  in  Fig.  2, 

has  been  found  to  approximate  the  aperture  region  A. 

We  associate  M  with  the  nth  edge.  As  Fig.  3  shows,  there  are  two 

triangles,  T^  and  T  ,  related  to  the  nth  edge  (assumed  not  on  the  boundary) 
n  n 

of  a  triangulated  area  modeling  the  aperture.  The  global  position  vector 

r  and  the  local  position  vectors  p^,  p  are  defined  as  shown.  The  plus 
—  “n  — n 

or  minus  designation  of  the  triangles  is  determined  by  the  choice  of  a 

positive  current  reference  direction  for  the  nth  edge,  which  is  assumed 

to  be  from  T^  to  T  .  Define 
n  n 

^  I 

■1  { 


Fig.  2.  Triangulation  Example. 


Fig,  3.  Expansion  function. 
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— T  P„  »  £  in 

2A- 

n 


0  *  elsewhere 


+  + 
where  H  is  the  length  of  the  edge  n  and  A  is  the  area  of  triangle  T  • 
n  n  n 

As  pointed  out,  ^  is  related  to  the  nth  edge,  which  is  not  on 

the  boundary  of  the  aperture.  Since  the  magnetic  current  must  not  have  a 

normal  component  on  the  boundary,  we  need  not  define  basis  functions  for 

any  boundary  edges.  (See  the  reason  in  the  following  section.) 

Using  the  basis  functions  above,  we  see  that  all  edges  of  t"*”  and  T 

n  n 

are  free  of  magnetic  line  charges.  For  the  common  edge  n,  the  normal  com¬ 
ponent  of  the  magnetic  current  is  constant  and  continuous  across  the  edge 


(see  Fig.  4),  shown  as  follows: 


^  a'*' 

w+  ,,  _n_  _ =  1 

n,  normal 

n 


M  ,  normal 
n 


I  A 
n  n 

2a'' 

n 


Hence,  V  may  be  interpreted  as  the  normal  component  of  the  magnetic 
n 

curr  i  .t  density  crossing  the  nth  edge.  For  the  other  conjoined  edges, 

M  has  no  component  normal  to  them,  and  hence  no  magnetic  line  charges 
“n 

exist  along  those  edges. 

With  basis  function  M  defined  as  above,  the  associated  magnetic 
— n 

surface  charge  density  is  of  the  form  of  pulse  doublets: 


m  =  -  V  •  M 
n  ju)  — r 


1  a  ”n) 


T  .  £  in  t; 


(8) 


Also,  It  can  be  proved  that  a  superposition  of  the  basis  functions 
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Y  =  4ja)  (W  •  F  +  'll  w  )dS 
mn  /  — m  — n  n  m 


4  {ju)£  [F  •  =f-+  F  (l!")  =f-]  +  i  [0  (r'"")  -  $  ] } 

m  — n  — m  2  — n  — m  2  m  n  — m  n  — m  '  ^ 


-  «  (J“!F  (4'")  •  1-  +  F  (4')  •  -F  ♦  <4')  - » <4*)) 


where 


^  •  G(k.  r!',  r’)dS' 

n  4Tr  “n  —  ~tr  — 


*  ,  c±^  -1 

$  (r  )  =  7— 2 — 
n  — m  4Trjajij 


V'  .  M  (r’)  G(k,  r'^",  r')dS' 
J  s  — n  ~  ~m  — 


To  evaluate  F  (r  )  and  $  (r  "),  we  proceed  face  by  face  for  the 
— n  “in  n  “m 

sake  of  efficiency.  Now, let  us  look  into  the  case  shown  in  Fig.  5,  with 

observation  triangle  T  and  source  triangle  T  .  The  number  of  basis 

P  q 

functions  for  T  is  less  than  or  equal  to  three. 

q 

p^  =  ±  (r '  -  r^) 


fP**  =  F  (r'^) 
“1  -q^  -p 


$5*^  =  0  (r*^) 

i  q^^  “p 


where  1  =  1,2,3.  Hence, 


i&JSl!;-.-'- 
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-1  Att 


,  M  (r')  G(k,  r')dS’ 

J  -q.  -  -p  - 


Att 


q 


i  =  1,2,3 


+1 

ATTjojy  ] 


G(k,  r^,  r’)ds’ 

\  -p  - 


i  =  1,2.3 


Now,  make  use  of  the  area  coordinate  [8]  for  triangle  T  .  (Check 

q 

Fig.  6.) 


i.e. 


where 


r '  =  ^rj^  +  nr^  +  Cr^ 

(x,y)  =  +  n(x2,  y2)  +  ^(x^,  y^) 

e  ■ 

n  =  A,/A 
2  q 

C  =  =  1  -  C  -  n 


We  transform  the  surface  integrals  for  and  into  double  Integrals  by 
the  following  formula 


1  1-n 

f(r’)dS’  -  2A^  I  J  f(Crj^  +  ^2  +  (l-C-n)r3)dCdn  (12) 

T-  0  0 

q 


r 
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(15) 


15 


Hence 


=  2  W  •  dS 
m  -m  — t 


where 


c+  c- 

-  •  =^  + 

m  t  2  ^  “TO  2 


=  (U  .H^  +  U  ,Hh  e 
-t  -m  9  ^ 


.,  1  c± 
3k  •  r 


e 


=  -  U  .k^ 
—  —  1 
r 


-  k^(x  sin  0^  cos  <!)^  +  ^  sin  6^  sin  ^  cos  6^) 


By  a  similar  formulation,  we  have 


(16) 


h“°  dS 


,vC+  c- 

.  P  P 

r„roO/  c+,  ~n  .  ..mo,  c-,  — n  , 

=  -  2iJ^  (I„  )  '“  +  H,  (r^  )  •  -2-1 


(17) 


where 


m  c± 

•"  1 1c  ’IT 

..ino  /  ci^  /tt 

(r  )  “  (U  Hq  +  U  H.)e 
— t  -n  — .m  9  —  m  4> 

9  ({) 

k"*  =  -  U  k"* 

—  —  m 

r 

*  -  k  (x  sin  6  cos  (p  +  ^  sin  0  sin  $  +  ^  cos  0  ) 


VII.  REPRESENTATIVE  QUANTITIES  TO  BE  CALCULATED 

We  first  calculate  the  equivalent  current  M  »  J  V  M  ,  then 

—  n  "n 

-1  ” 

the  equivalent  charge  density  m  =  —  '7  •  M,  the  far-field  the  inci¬ 
dent  power  V,  ,  and  the  power  transmitted  .  We  can  find  the 

inc  trans 

transmission  cross  section  patterns  t,  the  transmission  coefficient  T, 
and  the  transmission  area  TA.  They  are  listed  as  follows:  f'  ost  of  the 
derivations  are  in  [IJ.) 
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V  =  [Y®  + 


=  ^  [Y^"]-^  .  2 


(18) 


m  = 


+  ^  + 

I  7-1  ’ 

I  JwA 


0  ,  elsewhere 


H  =  ^  P™  [Y^®]-^ 

m  Srrr  ‘  ^ 


iwi 


87rr 


m 


2  i”°  I  [Y^^]-^  .  2 


-»i  Iqp 

-Icoe  g  m  jmo  ^lo 

4Trr  ®  ^  ‘ 

m 

T  *  2Trr^  •  ri|H  |^/ri|H^*’l^ 
in  ‘  m 


(19) 


Pine  °  S  cos  6^ 


^  ha 

P.ran=  =  Re(V[Y'‘®]  V 

trans 

1  ^  -*-i* 

=  Re(|  V  •  I-"  ) 

1  -V  -Vl* 

Re(V  I  ) 

ri|H^°|^  S  cos  0^ 


(20) 
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TA  =  T  •  S  cos  0^ 

=  Re(V  I^*)/2n|H^°|^  (21) 

VIII.  NUMERICAL  RESULTS  AND  DISCUSSION  FOR  THE  PROTOTYPE  CASES 

Using  the  previous  formulation  and  adopting  quite  a  few  subroutines 
from  [3],  we  have  developed  a  versatile  computer  program.  This  program  can 
solve  not  only  the  prototype  problems,  but  also  both  extensions  mentioned 
in  Section  I.  It  is  described  and  listed  in  Section  XIV.  Some  representa¬ 
tive  computations  for  the  prototype  cases  are  given  in  this  section.  To 
ensure  the  validity  of  our  formulation,  we  examined  several  special  examples 
which  are  available  in  the  literature.  As  we  will  see,  the  results  agree 
very  well. 

The  first  example  is  a  narrow  slot,  width  X/20  and  variable  length  L, 
lying  on  the  x-y  plane  with  axis  in  the  y  direction.  This  slot  aperture  is 
illuminated  by  a  normally  incident  plane  wave  with  unit  magnetic  field 
polarized  in  the  (()  direction.  Figure  7  shows  the  configuration.  As  shown 
in  Fig.  8,  it  is  triangulated  into  40  patches.  Figure  9  shows  the  trans¬ 
mission  cross  section  patterns  in  two  principal  planes,  i.e.  at  (()>  =  90“ , 

0  =  90  180“)  and  (*  =  270“,  0  =  180“  ->■  90“),  x .  at  (ij)  =  0“ ,  0  =  90“  ->  180“) 

<P 

and  (((>  =  180“,  0  =  180“  -*■  90“).  Figure  10  plots  the  equivalent  magnetic 
current  in  the  aperture  region.  Table  1  shows  the  physical  area  A,  trans¬ 
mission  coefficient  T,  transmission  area  TA  of  the  corresponding  cases. 

From  our  data,  we  found  good  agreement  with  earlier  results.  The  only  case 
of  significant  discrepancy  is  the  far-field  for  L  “  X/4,  and  it  is  believed 
that  our  result  is  correct. 
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Fig,  7. 


Slot  aperture  under  unit  H 


1 


normal 


inc idence. 


X 


1 


Fig.  8.  Tr iangulation  of  the  slot  aperture. 


Transmission  cross  section  for  slot  aperture  of  Fig.  7 
W  =  >720;  L  =  >/4,  >/2,  3>/4,  >  In  (a),  (b),  (c),  (d) . 
2 

+  :  computed  'f,7>  In  principal  plane. 

^  2 

computed  in  principal  plane, 

corresponding  results  from  [2]. 


Equivalent  magnetic  currents  for  slot  of  Fig.  7. 

M  is  in  y-direction. 

W  =  A/20;  L  =  A/4,  A/2,  3A/4,  A  in  (a),  (b) ,  (c),  (d) 
+:  computed  magnitude  of  |m/E  | . 

computed  phase  of  |m/^  |. 

0,/'.:  corresponding  results  from  [2]. 
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Table  1.  Transmission  coefiicient  for  slot  of  Fig.  7,  W  =  //20. 


L  =  X/4 

A 

T 

TA 

0.12500E  -  01  (X^) 

0.17950E  +  00 

0.22438E  -  02  (X^) 

A 

0.25000E  -  01  (X^) 

L  =  X/2 

T 

0.81829E  +  01 

TA 

0.20457E  +  00  (X^) 

A 

0.37500E  -  01  (X^) 

L  =  3X/4 

T 

0.21401E  +  01 

TA 

0.80254E  -  01  (X^) 

A 

0.50000E  -  01  (X^) 

L  =  X 

T 

0.15163E  +  01 

TA 

0.75815E  -  01  (X^) 

The  second  check  is  made  for  a  square  aperture  lying  in  the 
x-y  plane.  It  is  illuminated  by  a  normally  incident  plane  wave  with  unit 
magnetic  field  polarized  in  the  (j)  direction.  This  square  aperture  is 
triangulated  into  32  patches  as  shown  in  Fig.  11.  The  computed  far-field 
patterns  match  almost  perfectly  the  results  from  [2] ,  which  are  plotted 
in  Fig.  12,  For  near  field  quantities,  there  is  no  available  data  to  com¬ 
pare  with.  Nevertheless,  our  results  shown  in  Fig.  13,  seem  to  be 
reasonable. 

For  a  third  check,  we  consider  a  circular  aperture  which  is 
triangulated  into  24  patches  as  shown  in  Fig.  14.  There  are  two  things 


11*  Tr iangolation  of  the  square  aperture. 

Illuminated  by  a  normal  incident  plane 
wave  with  unit  magnetic  field  polarized 
in  (})-dlrection. 


A  =  0. 62500  E- 01 


T  =  021483  E+00 


TA=  0.13427 E-OI  (X*) 


A  =  0.25000  E+00(X*) 


T  =  0.15673  E+OI 


TA=  0.39183  E+00(X*) 


Fig.  12.  Transmission  characteristics  for  square  aperture 
in  Fig.  11 . 

All  the  notation  and  marks  are  similar  to  Fig.  9 
and  Table  1.  But  L  =  >/A,  )i/2  in  (a),  (b)  . 
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to  notice.  First,  to  compensate  for  the  loss  in  total  area,  we  9'  'id  put 
the  boundary  vertices  outside  the  cin  le  so  that  we  get  the  correct-  total 
aperture  area.  Second,  to  take  care  of  the  edge  effect,  we  need  a  higher 
patch  density  around  the  boundary  than  in  the  center.  Now  this  circular 
aperture  is  excited  by  an  obliquely  incident  plane  wave  with  either 
parallel  or  perpendicular  polarization  [3].  To  compare  with  the  data  available, 
we  redefine  the  transmission  coefficient,  denoted  as  TCHA,  Instead  of  the 
previous  T.  It  is  normalized  with  respect  to  the  incident  power  density 
at  normal  incidence  rather  than  the  actual  incident  power  density  at 
oblique  incidence.  As  can  be  seen  from  Fig.  15,  our  computed  data  agree 
well  with  the  previous  literature  [3].  The  slight  discrepancy  probably  is 
due  to  the  edge  effect  plus  the  difficulty  in  matching  the  exactly  circular 
boundary  with  straight  line  segments. 

So  far,  all  the  cases  we  have  tried  are  just  validity  checks. 
Obviously,  for  any  rectangular  aperture  (including  the  narrow  slot,  the 
square  aperture,  etc.),  triangular  patching  is  not  superior  to  the  rec¬ 
tangular  patching  in  [2].  Also,  for  any  annular  aperture  (including 
circular  aperture),  triangular  patching  is  not  superior  to  the  annular 
subsections  used  in  [3]. 

To  show  the  versatility  of  the  triangulation  method,  we  try  two 
other  shapes.  These  are  a  diamond-shaped  aperture  and  a  cross-shaped 
aperture.  Our  formulation  treats  them  without  difficulty  .  Figures  16 
and  17  demonstrate  some  simple  ways  to  triangulate  the  diamond  aperture 
and  cross  aperture  into  12  patches  and  20  patches. 
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IX.  EXTENSION  I:  HALF  SPACES  WITH  DIFFERENT  MEDIA 

Assume  two  half  spaces  with  different  media,  and  C^,  separated 
by  an  Infinite  conducting  plane  with  an  arbitrarily  shaped  aperture.  The 
method  of  solution  will  be  almost  the  same,  except  for  some  minor  dif¬ 
ferences. 

f  s 

First,  instead  of  Y  =4  <-W  ,  H  (M  )>,  the  elements  in  the 
mn  m  t  — n 

admittance  matrix  will  be 


Y^  +  Y*’  =  2<-W  ,  H^®(M  )  +  )> 

mn  mn  m  t  “n  t  — n 


cf 

=  2)1  F^r'^-)  •  =?-]  +  ^rcr"-)  - 


c+^ 


m  “H 


c+ 

+  2£  {jaj[F^(r^'^)  •  %-  +  F^(r^"'  ' 
m  ~~n  “m  z  — ti  — m 


c- 

m_ 

2 


+  F^r‘=")  •  (22) 


where 


a 

b  r 


pb  (,C±  ^  ^ 

— n  — m  UTT 


M  (r’)  •  G(k‘’  ,  If*^'  -  r’  |)ds' 


,.b  ,  c±,  -1 

V  (t  )  =  — I — 

n  — m  ATTjuM 


V'  •  M  (r’)  G{k^  ,  Ir'^'  -  r’|)ds' 
s  “n  —  "Tn  — 


Hence 


e’’  £ 


pPq  .  + _ i  (j.  tP'I  +  r  iP'l  +  r  iP**  -  r  iP**) 

^  i  Att  ^^1  -2  S  ^  ^3  \  -i  ^  ^ 


where 


a 

b 


'a 

b 


/e 
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Second,  use  =  i/e  for  the  excitation  vector  I^,  and  kT  =  /e.  k™ 
—a  a  —  b  — 

for  the  measurement  vector  Y".  Then  all  the  remaining  calculations  are 
exactly  the  same  as  in  the  prototype  problem. 

X.  NUMERICAL  RESULTS  AND  DISCUSSION  FOR  EXTENSION  I 

Some  arbitrary  but  interesting  combinations  of  different  e  and  e. 

a  b 

are  used  as  examples.  Tables  4,  5,  and  6  give  the  peak  values  of  the 
dominant  component  of  M,  the  transmission  coefficients,  the  transmission 
areas,  and  the  maxima  of  the  transmission  cross  section  patterns  in  the 
principal  planes. 

The  sources  are  normally  Incident  plane  waves,  with  the  magnetic  field 
polarized  along  the  largest  dimension  of  the  aperture,  for  the  diamond-shape, 
cross-shape,  circular,  etc.  All  the  largest  dimensions  are  chosen  to  be  a 
quarter  wavelength,  and  the  relative  dielectric  constants  are  chosen  to  be 
combinations  of  1  and  4.  The  results  show  an  interesting  phenomenon;  i.e. 
after  changing  the  dielectric  constant  on  one  side,  the  aperture  appears 
to  be  resonant  with  respect  to  that  half  space.  Hence  both  the  equivalent 
current  and  the  transmission  characteristics  have  significant  increases. 
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Table  A.  Characteristic  quantities  of  a  diamond  aperture  with  respect 

to  some  combinations  of  and  e, 

^  b 

Aperture  area  =  0,0036  m^,  L  =  0.25X  ,  W  =  L/2  where  X  is  the 

o  o 

free  space  wavelength  of  the  incident  plane  wave.  (X^  =  0.48m) . 


K-  s' 

(1,  1) 

(1,  4) 

(4,  1) 

(4,  4) 

0.58000E  00 

0.91000E  00 

0.72000E  00 

0.11450E  01 

0.14300E  01 

0.22890E  01 

0.18000E  01 

0.28680E  01 

T 

0.23621E  -01 

0.15036E  00 

0.12028E  01 

0.33212E  01 

TA<‘'> 

0.36908E  -03 

0.23493E  -02 

0.18794E  -01 

0.51892E  -01 

0.57259E  -03 

0.56648  -01 

0.88513  -03 

0.89017E  -01 

0.57259E  -03 

0.14162E  -01 

0.22128E  -03 

0.89017E  -01 

O.OOOOOE  00 

O.OOOOOE  00 

O.OOOOOE  00 

O.OOOOOE  00 

0.55987E  -03 

0.12919E  -01 

0.21629E  -03 

0.81170E  -01 
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Table  5.  Characteristic  quantities  of  a  cross  aperture  with  respect 

to  some  combinations  of  e  and  e,  . 

a  b 

2 

Aperture  area  =  0.008  m  ,  L  =  0.25A  ,  W  =  L/3  where  X  is  the 

o  o 

free  space  wavelength,  (X^  =  0.48  m) . 


(1.  1) 

(1,  4) 

(4;  1) 

(4,  4) 

|m/e| 

'  ‘max 

0.78000E  00 

0.11280E  01 

0.22560E  01 

0.19000E  01 

T 

0.89157E  -01 

0.79526E  00 

0.63621E  01 

0.84532E  01 

0.30957E  -02 

0.27613E  -01 

0.22091E  +00 

0.29351E  +00 

0.49090E  -02 

0.71640E  00 

0.11194E  -01 

0.54156E  00 

0.49090E  -02 

0.17910E  00 

0.27984E  -02 

0.54156E  00 

O.OOOOOE  00 

O.OOOOOE  00 

O.OOOOOE  00 

O.OOOOOE  00 

0.46000E  -02 

0.14100E  00 

0.26300E  -02 

0.43650E  00 
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Table  6.  Characteristic  quantities  of  a  circular  aperture  with  respect 

to  some  combinations  of  e  and  e.  . 

a  b 

2 

Aperture  area  =  0.19635  m  ,  R  =  0.25X  ,  X  =  Im  (free  space) 

o  o 

2  ~ 

(Triangulized  aperture  area  =  0.17678  m  ,  X  =  0.94886  m.) 


(1.  1) 

(1,  4) 

(4,  1) 

(4,  4) 

|m/e| 

0.18860E 

01 

0.11212E 

01 

0.22420E 

01 

0.15430E 

01 

0.44320E 

01 

0.16350E 

01 

0.32700E 

01 

0.19740E 

01 

T 

0.15820E 

01 

0.89937E 

01 

0.71949E 

01 

0.46297E 

01 

0.27967E 

00 

0.15899E 

01 

0.12719E 

01 

0.81844E 

00 

(T^/X^) 

0.72745E 

00 

0.12550E 

02 

0.19610E 

00 

0.52419E 

01 

0  max 

irjxh 

0.72745E 

00 

0.31375E 

01 

0.49024E 

-01 

0.52419E 

01 

<j)  max 

O.OOOOOE 

00 

O.OOOOOE 

00 

O.OOOOOE 

00 

O.OOOOOE 

00 

0.32401E 

00 

0.68000E 

-01 

0.24698E 

-01 

0.14000E 

00 
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XI.  EXTENSION  II;  LOSSY  DIELECTRIC  WINDOW 

Assume  now  that  the  arbitrarily-shaped  aperture  is  covered  with  a 
non-magnetic  material  sheet  (u  =  which  is  very  thin.  Then,  we  can 
treat  this  aperture  as  a  continuously  loaded  case.  The  addition  of  one 
more  term  to  the  admittance  matrix  of  the  prototype  problem  will  suffice 
to  give  the  solution  in  a  straightforward  way.  Both  theoretical  and 
numerical  derivations  will  be  given  in  moderate  detail. 

Since  the  aperture  region  is  covered  with  a  lossy  dielectric 
(O,  c)  there  will  be  current  [10]  consisting  of  conduction  current  and 
polarization  current.  This  gives  a  total  increase  of  volume  current 
density 

A  J  =  J  =  oE  +  ^ 

—  — m  —  dt 

==a  E  +  jo)  (£  -  e^)E  (23) 

Since  this  window  is  assumed  to  be  very  thin,  the  electric  fields  on 
both  sides  are  still  continuous.  Hence  equivalent  currents  can  still 
be  M  in  region  a  and  -M  in  region  b.  The  current  has  no  normal  component 
and  is  tangential  to  the  window.  The  Increase  in  surface  current  density 
at  the  aperture  region  should  be 

AJ  =  d  J  d  <  A./10 

— s  — m 

where  d  is  the  thickness  of  the  window. 

Now,  Instead  of  as  in  the  prototype  problem  the  boundary 

condition  at  the  aperture  region  is 
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n  X  (h'’  -  H®)  =  AJ 
-  '  -s 

=  [a  +  j(j)  (e  -  e  )]dE 
o  — 

=  [a  +  j(jjA£)d  •  (-  n  X  M) 


where 

=  (a  +  ja)Ae)d  (24) 

In  general,  for  a  good  dielectric,  the  conduction  current  is  much 
smaller  than  the  polarization  current.  Then  =  jwAed,  i.e.  purely 
susceptive. 

Equation  (24)  can  be  rewritten  as  follows; 

]i^(-M)  -  hJ(m)  -  H J  =  ~  y^H 

-  -  h!(M)  +  y  M  = 

— t  —  — t  —  y—  — t 

’  a  h  y 

[Y  +  Y  +  Y  ]V  =  I  (25) 

where  j 

IV  1  -  c<w.. 

If  the  window  is  isotropic  and  homogeneous,  then 

tv‘l  -  V'm  «« 

and  ^  “  [Y^  +  Y^  +  Y^]  ^I^  can  be  solved  by  analogy  to  the  prototype 

S, 

problem  as  soon  as  we  have  Y  computed. 
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XII.  NUMERICAL  RESULTS  AND  DISCUSSION  FOR  EXTENSION  II 

As  a  first  example,  let  us  consider  some  good  dielectric  sheets 
-4 

(low  loss:  0  =  10  )  with  different  dielectric  constants  being  900, 

81,  and  4,  and  variable  thickness  d  ranging  from  zero  (uncovered  aper¬ 
ture)  to  some  reasonable  limits  (e.g.  around  O.lA,,  where  A,  =  A  //i~) . 

d  dor 

We  use  these  materials  for  diamond-shaped  windows  (major  axis  0.45A  , 

o 

minor  axis  O.IX^).  Each  is  illuminated  by  a  normally  incident  plane  wave. 

On  first  thought,  we  might  expect  the  transmission  coefficient  to 
become  smaller  after  we  cover  the  aperture  with  a  lossy  dielectric  window. 
However,  in  the  limiting  case  (magnetic  dipole  mode  for  a  small  aperture), 
the  susceptance  of  a  small  aperture  is  inductive  [17].  Now  since  the 
dielectric  sheet  we  use  is  essentially  a  distributed  capacitive  loading, 
we  might  expect  some  kind  of  "resonance-like"  behavior  to  occur.  Our 
results  support  this  expectation. 

Figures  19,  20,  and  21  show  the  transmission  coefficients  vs.  thick¬ 
ness  for  the  diamond-shaped  windows.  The  incident  wavelength  A^  is  0.2m. 

We  see  from  these  figures  that  when  we  start  increasing  d,  the  transmission 

_2 

coefficients  drop  from  their  original  value  (0.504)  to  almost  zero  (-10  ). 

But,  then,  instead  of  becoming  exactly  zero,  there  is  a  jump  in  each  case. 
This  resonance-like  behavior  can  be  explained  as  the  result  of  the  better 


match  between  the  two  half  spaces  provided  by  the  dielectric  sheet  of 

proper  thickness.  Even  though  this  resonance  occurs  at  different  d's 

for  different  materials,  it  always  reaches  the  same  maximum  value  (8.728), 

i.e.,  the  T  is  Independent  of  the  material  (e  ).  Also,  this  T 

’  max  r  max 


is  less  than  the  optimum  value  for  small  apertures  (T^p^  ~  “  10.61, 


0  2.0  4.0  6.0  8.0  DO  12.0  14.0  160  180  20.0  UlO'®) 

d 

Fig.  19.  Transmission  coefficient  of  a  diamond-shaped  window,  e  =  900. 

,  -4 

=  0.2;  o  =  10  . 


0  2.0  40  6.0  8.0  10.0  120  14.0  ISO  18.0  20.0(xl0-<) 

d 

Fig.  20.  Transmission  coefficient  of  a  diamond-shaped  window,  e  =81. 

A  =  0.2;  a  =  lo“^. 
o 
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at  resonance).  This  is  due  to  the  fact  that  the  capacitive  coupling  through 

the  dielectric  sheet  can  never  reach  exactly  perfect  matching,  since  we  are 

only  modifying  five  elements  of  each  row  in  the  admittance  matrix. 

-4 

As  a  second  example,  we  choose  o  =  10  ,  =  900  for  a  rec¬ 

tangular  window  (actual  slot  0.45^^  by  0.05X^)  illuminated  by  a  normally 
incident  plane  wave.  Figure  22  shows  the  result.  Since  the  slot  is  already 

in  resonance,  its  transmission  coefficient  has  the  optimum  value  (T  = 

max . 

10.8  =  T  )  before  we  cover  it.  Then  T  monotonically  decreases  as  we  In- 
opt . 

crease  d.  There  is  a  small  rise  near  d  =  0.5X,,  but  that  thickness  is 

d 

probably  outside  the  range  of  our  theory.  We  can  say  that  the  window  pro¬ 
vides  a  shielding  effect  in  this  case. 

Finally,  let  us  check  the  effect  of  high  loss  dielectric  materials. 

With  the  same  rectangular  window  and  the  same  incidence  as  in  the  previous 

example,  we  choose  four  different  cases:{a,  c^)  =  (5000,  900),  (5000,  4), 

4  4 

(10  ,  900),  (10  ,  4).  The  results  are  shown  in  Figs.  23  and  24.  These  curves 

are  interesting  yet  difficult  to  interpret,  especially  for  the  extremely  high 

and  narrow  peak  right  before  the  dropping  to  zero.  Regardless,  they  show 

very  good  shielding  effect;  l.e.  the  T  already  drops  to  zero  at  around 
-4  -4 

d  *  2 X  10  ,10  .  Actually,  since  the  skin  depths  here  (approximately 

1.84  X  10  ^  and  1.30x  10  ^)  are  much  smaller  than  one  tenth  of  X^  (6.67x10  ^ 

_2 

and  10  ),  they  play  a  more  Important  role.  As  we  can  see,  the  unexplanable 

peaks  occur  near  those  thicknesses.  Tt  may  be  that  our  formulation  doesn't 
work  when  the  sheet  thickness  become  comparable  to  its  skin  depth.  Never¬ 
theless,  with  high  conductivities,  changing  dielectric  constant  doesn’t  change 
the  curve  much;  while  with  the  same  dielectric  constant,  doubling  the  conductivity 
gives  the  same  effect  with  only  half  of  the  thickness  required.  This  is  reason¬ 
able  because  there  the  factor  a  +Si»hz  is  approximately  a  +j(e^-l)/12,  which 
obviously  has  a  dominant  real  part. 


XIII.  CONCLUDING  REMARKS  | 

I 

To  obtain  the  transmission  characteristics  of  arbitrarily-shaped 

i 

apertures,  the  generalized  network  formulation  for  aperture  problems  was 

1 

utilized.  Triangular  patching  and  local  position  vectors  were  used  as 
bases  for  arbitrary  2-dimensional  shapes.  Extensions,  such  as  different 
media  in  half  spaces  and  lossy  dielectric  windows,  were  derived.  Programs 

for  calculating  both  far-field  and  near-field  quantities  were  also  developed.  ■ 

The  apertures  considered  in  this  report  were  for  coupling  between 
two  half  spaces.  As  further  work,  we  could  try  some  of  the  following: 
arbitrarily-shaped  apertures  backed  by  an  inf initely-long  wire,  and  arbitrarily- 

/ 

shaped  aperture  providing  coupling  between  various  combinations  of  half-spaces-,  ■ 
waveguides,  and  cavities. 

We  should  point  out  that  there  are  systematic  schemes  for  generating 
the  triangular  patching  model  [14,  15].  But  for  the  sake  of  simplicity,  it 
may  still  be  preferable  for  the  user  to  input  all  the  nodes  and  meshes  himself. 

One  more  thing,  magnetic  and  electric  polarizabilities  of  small  aper¬ 
tures  have  been  solved  with  scalar  bases  and  quadrilateral  and  triangular 
patches  [16].  It  would  be  interesting  to  find  the  polarizabilities  for 
arbitrarily-shaped  small  apertures  by  using  the  vector  bases. 


XIV.  COMPUTER  PROGRAMS 

This  complete  program  treats  the  most  general  case;  arbitrarily-shaped 
apertures  on  an  Infinite  conducting  plane,  half  spaces  on  both  sides  can 
have  different  media,  the  aperture  can  be  covered  with  a  lossy  dielectric 
sheet.  A  thorough  listing  of  the  complete  program  will  be  given  after  a 


short  description. 


A9 


In  addition  to  the  main  program,  there  are  18  subroutines  and  1 
function  subprogram  in  this  set.  They  are  the  following: 

MAIN 

SOLTN 

INDATA 

GEOM 

AJUNC 

CURDIR 

BODPAR 

YMATRX 

YWINDO 

MAGCHA  / 

TRANS 

MEASUR 

SCAINT 

VECINT 

LININT 

INTGRL 

CA 

EXPRN 

CSMINV 

DTRMNT 

Brief  descriptions  of  MAIN,  SOLTN,  AJUNC,  YMATRX,  YWINDO,  MAGCHA, 

TRANS,  MEASUR  will  be  given;  all  the  other  subroutines  were  adopted  from  Rao's 
work  and  modified  for  2-dimensional  geometry. 
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MAIN  program  reads  in  the  number  of  nodes  NNODES ,  the  number  of 
edges  HEDGES,  and  dielectric  constants  DIA  and  DIB.  It  calls  subroutine 
INDATA  to  read  in  all  the  nodes  DATNOD  and  meshes  NCONN.  By  calling 
subroutines  GEOM,  AJUNC,  CURDIR,  BODPAR  it  obtains  the  triangular  patching 
model  of  the  aperture.  Then  subroutine  SOLTN  is  called  to  find  the  trans¬ 
mission  characteristics. 

SOLTN  subroutine  reads  in  properties  of  the  possible  window;  con¬ 
ductivity  SIGMA,  difference  between  the  dielectric  constant  and  1  (free 
space)  DELDIE,  thickness  THICK;  window  with  zero  thickness  means  aperture 
without  covering.  It  also  reads  in  the  number  of  incident  fields  to  be 
treated  NFIELD,  incident  angles  THETA  and  PHI,  incident  polarizations 
ETHETA  and  EPHI,  and  a  flag  IRCS  for  calculating  the  aperture  transmission 
cross  section.  It  calls  subroutine  YMATRX  and  YWINDO  to  get  the  admittance 
matrix  CY  and  the  excitation  vector  Cl,  inverts  CY  by  calling  CSMINV,  then 
obtains  the  magnetic  current  coefficients  CV  and  calculates  the  transmission 
area  TS,  transmission  coefficients  T  and  TCHA  (the  latter  is  normalized  with 
respect  to  normal  incident  power  density).  As  an  option,  subroutine  MAGCHA 
can  be  called  to  give  the  corresponding  magnetic  charge  distribution.  Finally, 
IRCS  serves  as  an  indicator  to  show  whether  subroutine  TRANS  should  be  called 
and  which  polarizability  of  the  transmitted  field  is  to  be  considered; 

IRCS  =  0  means  the  TCS  pattern  is  not  desired,  1  or  2  means  theta-  or  phi- 
polarized  TCS  pattern  is  to  be  computed. 

AJUNC  subroutine  finds  the  corresponding  pair  of  triangles  for  each 
non-boundary  edge,  NJUNC.  It  serves  as  a  preparation  for  YWINDO. 

YMATRX  subroutine  calculates  the  summation  of  admittance  matrices, 
for  both  half  spaces, CY,  the  excitation  current  vector  Cl  in  half  space  A. 
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It  evaluates  the  necessary  integrals  in  the  triangular  area  by  calling 
SCAINT  and  VECINT. 

YWINDO  subroutine  is  called  only  when  THICK  is  not  zero.  It  com¬ 
putes  the  load  admittance  and  adds  it  to  CY. 

MAGCHA  calculates  the  equivalent  magnetic  charge  densities.  It 
can  be  dropped  out  without  hurting  the  completeness  of  the  general-purpose 
program  set . 

TRANS  subroutine  is  called  to  find  the  theta-  or  phi-polarized  TCS 
pattern  when  IRCS  is  1  or  2.  It  reads  in  the  initial  angles  THETAl  and 
PHIl,  final  angles  THETA2  and  PHI2,  numbers  of  subangles  NTHETA  and  NPHI 
to  form  the  mesh  nodes  for  desired  TCS.  It  calls  MEASUR  to  obtain  the 
measurement  vector  CIM  in  half  space  B. 

MEASUR  subroutine  computes  the  measurement  current  vector  CIM  in  a 
similar  way  as  YMATRX  computes  Cl. 
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C  THIS  PROGRAM  CALCULAT6S  THE  MAGNETIC  CURRENT  DISTRIBUTION 
C  FOR  A  ARBITRARY-SHAPE  APERTURE  CN  A  PERFECT  CCNOUCTING  PLANE 
C  EXCITED  BY  A  PLANE  WAVE  OF  DESIRED  AMPLITUDE  E-FIELO. 

C  THE  ARBITRARY  APERTURE  IS  CESCPIBEO  BY 

C  A  SET  CF  NODES  AND  THEIR  COORDINATES. AND  EDGES  CONNECTING 
C  THESE  NOOES.THIS  PRCGRAM  MANIPULATES  THE  GIVEN  DATA  TO 
C  DETAIN  TRIANGULAR  PATCHES  AND  CURRENT  DENSITY  IS  CALCULATED 
C  AT  THE  CENTER  CF  EACH  EOGE.AS  IT  IS  THIS  PRCGRAM  CAN 
C  HANDLE  100  X  100  UNKNOWNS.  BUT  BY  MODIFYING  THE  DIMENSION 
C  STATEMENTS  CNE  CAN  USE  STILL  LARGER  SI2E  MATRICES. 

C  THE  MAIN  PRCGRAM  READS  THE  TCTAL  NUMBER  OF 

C  NODES  AND  TOTAL  NUMBER  OF  EOGES.IT  ALSO  CALLS  FIVE  SUBROUTINES. 
IMPLICIT  COMPLEX  IC) 

COMPLEX  CVI200).CI(200I.CY(SO.S0I.CIM(200> 

DIMENSION  OATNODI 1S0.3I 

INTEGER  NCCNNOOO.S).  ITRAK  (300  I.NBOUNDI  ISO.  4  )  .  I M 1  N(  300 1 

INTEGER  NJtNC(3GQ.3t 

COMMON/IF/IFACE 

CCMMCN/DIELEC/CIA.CIB 

C  READ  THE  TOTAL  NUMEEP  CF  NOCES (NNCOES )  AND  ECGESINEDGES I. 

READd.AS)  NNOOES.NEOGES 
45  FCRMAT(213I 

CALL  INOATA(OATNCC.NCONN.NNCOES.NECGES) 

CALL  GECM(NCONN.NEOUNO. ITRAK . IMIN .NE06ES . NNCOES . NSE .NAPTRE . NHANOL I 
C  HERE  WE  CALULATE  THE  TOTAL  NUMBER  OF  FACES! NFACES ) . 

C  WE  ALSC  CALCULATE  THE  NUMEER  OF  UNKNOWNSINUNKNS >• 

C  THE  NUMBER  OF  UNKNOWNS  ARE  OBTAINED  BY  SUBSTRACTING  THE 
C  NUMBER  CF  SURFACE  ECGESINSEI  FROM  TOTAL  NUMBER  OF  EDGES. 

NFACES=IFACE 

NUNKNS«NCOGES-NSE 

CALL  AJUNCINJUNC.NBQUNO.NECGES.NFACESI 
CALL  CURDIRINCCNN.NBOUNO.NFACES.NEOGES. IMIN.NSEI 
SO  CALL  BOOPARIOATNCO.NCONN.NBCUNO.NNOCES.NEOGES.NFACES. 

TNUNKNS. NSE. NAPTRE. NHANOL I 
REAOd.SSI  OlA.OlE 
SS  F0RMAT(2Fe.3) 

WR1TE(3.«01  OIA.DIB 

60  FORMAT! • I • ./// IX . • 0  I A= • .F8.3.//. tX. *0 IB= • .FE.3 I 

CALL  SCLTNICY. CV. Cl .NUNKNS.DATNOO.NCCNN.NCCUNO. HEDGES. 
SNFACES.NNQOES. IMIN. ITRAK. C IM.NJUNCI 
STCP 
END 


SUBROUT  INE  SOLTKCCYtCV  aCl •KUKKRSvO^TKCO* ^CC^N•Ne0U^0•NE0GeS• 
SNFACESvNNOCES* IMIN. I TR AK •€ IR .K JUMC I 
C  IN  THIS  SUBRCUTINE  .THE  MATRIX  ECUATION  VVsf  IS  SOLVED. 

C  V-NATRIX  AND  I-MATRIX  ARE  CBTAINEC  BV  CALLING  THE 
C  SUBROUTINE  VMATRX.  USING  SUBROUTINE  CSMINV.  HE  INVERT 
C  THE  V-MATRIX  AND  THEN  MULTIPLY  EY  CURRENT  VECTOR  TO  GET  THE 
C  VOLTAGE  COEFFICIEN1S.THIS  SUBRCUTINE  ALSC  C*LLS  T«0  OTHER 
C  SUBROUTINES.  NAMELY  MAGCHA  AND  TRANCS.  TO  CCMPUTE  MAGNETIC  CHARGE 
C  OISTRIEUTICN  ANC  TRANS  CROSS  SECTION  RESPECTIVELY. 

IMPLICIT  CCMPLEA  ICI 
REAL  CABS.COS 

COMPLEX  CY (NUNKNS.NUNKNS).CV4NUNKNS>.CICNUNKNS>.CIM(NUNKNSI 
COMPLEX  HTHETA.HPHI.ETHETA.EPHI.HMT.HMP 
DIMENSION  OATNOOINNOOES.ai 

INTEGER  NCCNN(NECGES.3).NB0UNC( ISO. 4 I. I TRAKI HEDGES I • I MINI  HEDGES > 

INTEGER  NJLNCIHEOGES.SI 

COMMON/PAR AM/THETA .PH  I . IF lELO 

CCMMCN/FIELO/ETHETA.EPHl.ALAMCA 

CCMMCN/MOOl/ARE AT .R AV . I KMAX • IMM IN .RM AX • RM IN 

COMMON/MODE/AVAREA. JMAX.ARMAX. JMIN.ARMIN.KMIN.RATIO 

COMMCN/FHINCI/FHINC 

COMMOK/KKK/AK.PI 

C0MM0N/0IELEC/0IA.0I8 

COMMCN/PCLARM/HMT.HMP 

COMMCN/FRE/ACMEGA 

C  SPECIFY  THE  VAVE  LENGTH  lALAMCAl  IN  METERS. 

PI=3.14ISS265 

ALAMOAsO.a 

OlEaSQRTf OlAI 

ALAMCAsALAMOA/0  IE 

AK«2«0*PI/ALAMCA 

VEL«3.0EHOE/01E 

FREQs!  VEL/ALAMOA  MI  .OE«0« 

ACME6AS2.04PI41.0EH064FRE0 
AMUs4.0*PI«I.CE-a7 
AIMPst20.0«Pl/0IE 
EPSLCNsl .0/<VEL4*2*AMUI 
AREAT=AREAT/f ALAM0A442) 

RAVsRAV/ALAMOA 
RMAXsRMAX/ALAMCA 
RMINsPMIN/ALAMCA 
AVAREAsAVAREA/I ALAM0A442) 

ARMAXsARMAXy<ALAMCA«*2l 
ARMINsARMlN/I ALAMCA442 I 
FACESMsl.O/AVAREA 
IIRITE<3.20SI 

205  FCRMAT(*t<//////////25X.*MC0ELING  PARAMETER  LIST  R  REGION  A*///) 
•RITEI3.206)  AREA! 

206  FORMATC/IOX.'SURFACE  AREA  CF  THE  APERTURE*  ■ • tE13.S. 1 X. • SO. 

S  MAVE  LENGTHS* ) 

HRITEf3.20GI  RAV. IKMAX.RMAX.IKMIN.RMIN 

209  FORMAT|/IOX.*AVERA6E  EDGE  LENGTH**. IE  13. 5. IX. • HA VE  LENGTHS* • 
•//tOX. •MAXIMUM  EDGE  LENCTHIEOGE  NO.*. 13. IX. • >*•• IEI3.5. IX. 

S*«AVE  LENGTHS* //tCX.*MINIMUM  EDGE  LENGTHIEOGE  N0.*.I3. 

SIX.* 1«*. IE  13.5. IX.* WAVE  LENGTHS* I 

MR ITE <3.210 >  AVAREA. JMAX.ARMAX. JMIN.ARMIN 

210  FORMATI/IOX. 'AVERAGE  FACE  AREA  *• • IEI3.S. tX. *SO.HAVE  LENGTHS*. 
S//10X. •MAXIMUM  FACE  AREA  CFACE  NO. • .1 3. 1 X. * }** . IE I3.S. I X. 
S*SC.«AVE  lengths*. //10X.*MIN|MUM  FACE  AREA  IFACE  N0.*.I3. 
•IX.*I**.1C13.S.IX.*S0.MAVE  LENGTHS* ) 
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MRITEO.at  tl  K*lIN«RAriO«FACESM 

211  FQFMATI/10}i*«»«IMRUM  FACE  FEICHT  TC  EASE  RATIO  ( 

SFACE  NO.'.ia.lX**  )=*.1E13.S«//10X«*AVERAGE  NUMBER  OF  FACES  PER 
*  SCUARE  NAVE  LENCTh^* . 1 E 13 .3  I 

I|R1TE(3.212> 

212  FORMATC«1«.//20X.«ELECTR1CAL  PARAMETERS*) 

WRITE! 3*2 13)  FREO« AOMEGA.AK.ALAMDA 

213  FORMAT! /10X**FREQUENCYs*«IC13.S»lX«*MF2*.//« lOX, •ANGULAR 
«FREQUENCVs«,lE13.5a lX.*RA01ANSySEC*»//»10X«*WAVE  NUM8ER«*« IEI3*S» 
SI X,*! 1/METERS )*«y/. 10X» •WAVE  LENGTFs* • IE 13. 5. 1 X. 'METERS* > 

tlR|TE!3.2tA)  ERSLCN<AMU<VEL<AIWP 
21 A  F0RMAT!/10X<*PERM1TTI VITV«* <1E13.5<1X< •FAR AOS/METER* < 
S//<10X<*PERNEAeiLtTV3*< 1EI3.S< 1X< ■NENRVS/METER*  < 

S//<10X<*VELCCITY  CF  L IGNTa • • iei3.S< 1 X< • METERS/SECONO* < 

S//<10X< •INTRINSIC  IMPE0ANCEa*<lE13<S<lX<*CNMS* ) 

C  INITIALIZE  THE  V-MATRIX 
00  10  Ia|<NUNXNS 
00  10  Jsl<NUKKNS 
CV! I< J)aCMPLX! 0.0<0.0> 

10  CCNTINUE 

C  READ  THE  NUMBER  OF  INCIDENT  FIELDS  FOR  WHICH  THE  MAGNETIC  CURRENT 
C  OISTRIEUTION  NEEDS  TO  BE  COMPUTED.  IT  IS  NECESSARY  TO 
C  SPECIFY  THE  INCICENT  FIELC  PARAMETERS  ON  SEPARATE  CARDS 
C  SO  THAT  THE  PRCGRAM  EXECUTES  CNE  SET  OF  PARAMETERS  AT  A  TIME. 

READ! 1<699)SI6MA<0EL0IE<THICK 
699  FCRMAT!E10<3.F10.S.Eia.3 ) 

REA0<1< ISINFIELC 

15  FORMAT! 13) 

OC  499  lJal<NFI£LC 
IFIELOal J 

C  INITIALIZE  THE  CURRENT  ANO  VOLTAGE  VECTORS. 

OC  498  fal.NUNKNS 
CV(I)aCNPLX!0.0<0«0) 

C  I  !  I  )  aC  MPL  X  !  0 . 0  <  0  •  0  ) 

498  CONTINUE 

C  READ  THE  INCICENT  FIELD  PAR AWETERS.THET A  ANC  PHI  REPRESENT 
C  THE  USUAL  SPHERICAL  CCCROINATE  ANGLES  WHICH  DETERMINE  THE 
C  DIRECTION  OF  PROPAGATION  OF  THE  PLANEWAVE.  ETHETA  AND  EPHI 
C  REPRESENT  THE  AMPLITUGE  OF  THE  INCIDENT  PLANEWAVE. 

C  THE  VARIABLE  IRCS  IS  REFERRED  TO  THE  CCMPUTATICN 

C  OF  TRANS  CROSS  SECTION.  IF  IRCSaQ.  THEN  RADAR  CROSS  SECTION  IS 
C  NOT  CCMFUTEO. 

READ! 1<I6)THETA<PHI <ETHETA<EPHI<IRCS 

16  FORMAT!6E10.3< 13) 

WRITE!3<17)  THETA<PHI 

IT  FORMATI//5X<"ANGLE  CF  INCIDENCE* <//< 10X<*TMETAa« , IE13<5< 1X< •DEGREE 
SS*<//10X<*PHIa*. IE 13.S< |X< •DECREES* ) 

HPHIaETHETA/CMPLXIA IMP.O.OI 
HTHETAaEPHI/CMPLXf-AIMP<O.OI 

FHINCaCABS!CMPLX!CA8S!HPHI )<CA8S!HTHETA) > ) 

WRITE !3< 18)  ETHETA, EPHI <HTHETA,HPH I 

18  FORMATI//5X<«PCLARIZATICN»<//<IOX<»B-TMETAa  f*<2EI3.S< 

••)  VOLTS/METER •< 

•//<IOX<*E-PHIa  !*<2EI3.S<*1  VOLTS/METER • < 

*//,10X<*H-THETA*! • <2E13<S<* I  AMPS/MBTER», 

*//<10X< 'H-PHIal •<2E13.S<* )  AMPS/METER*  > 

THETA«THETA4P 1/180.0 
PH|aPH|4PI/ie0.0 

C  CALL  THE  YMATRX  SUROUTINE  TO  FILL  THE  Y-MATRIX 
C  ANO  THE  CURRENT  VECTOR. 
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CALL  YI«ATRX<CV.CATNQD»NCONK«KeOUNC*N»tOOES«NEOGES* 

SNFACES.KUNKNS. ITDAK.CI) 

ZEfiC=.000E«00 

IF  <JHICK,6C,ZePCi  6C  TC  ««« 

CALL  VH  INOCfCV.CATKao.NCQNN.NeOUNO.NNOOES.NEOGES.NFACES.NUNKNS. 

SI  YFAK.KJUNCvSIGI^AcOELOlE.TMCK) 

999  IFdJ.NE.ll  GC  TC  19 
C  OBTAIN  THE  INVERSE  OF  THE  V-MATRIX. 

CALL  CSMINVICV.NUNKNS.NUNKFS.COTRM.ACCNC. lER) 

c  multiply  the  v-ikverse  matrix  rith  current  cclumn  vector 

C  TO  OBTAIN  VOLTAGE  CCEFF lEC lENTS. 

19  CC  20  Isl.NUNKNS 
DC  20  Jsl.NUNKNS 

evil )-CV( ll^CVI 1. J lAClIJ) 

20  CCNTINUE 

C  MRITE  THE  MAGNETIC  CURRENT  DENSITY  TABLE 

EABSsS0RT(CABS<EThETAI«*2»CAeS(EPHI 1**2) 

*R1TE<3«22) 

22  FORMATI* 1* •///2SX.*SURFACE  MAGNETIC  CURRENTS* •/// ) 

»RITE(3.23 1 

23  FCRMATCIX. *£006  NUMBER* •2SX. •MAGNET IC  CURRENT  DENSITY  CHEBS/M-Sl*) 
ARITE<3«24) 

24  FORMAT(/20X«*REAL*«EX«* IMAGINARY**6X«*MAGNITU0E*  aBXs'PHASE* 
*«ax.*|M/Ei  RAT1C*«6X,*PFASEC* ) 

NSEaNEOGES-NUNXAS 

KlsQ 

CC  SO  Isl.NEOGES 
IFINSE.EG.O)  GC  TC  3l 
OO  30  J^l.NSE 
IF< I.EC.fTRAKlJl)  GO  TO  4S 

30  CONTINUE 

31  K13K14I 
RA1«REAL(CV<K1 1 ) 

RA2«A1MAG(CV<K1  1) 

RA3<CA6S<CV<K1 1  1 
RA4<ATAN2<RA2<RA1 1 
RAS<RA3/EA6S 
RAe3RA4*ieC«0/Pl 

«R1TE<3« 101 1  I vRA I, RA2. RA3<RA4<RAS<RAB 
101  F0RMAT<2X«14.8X<1E13.S<2X, 1E13<S<2X< 1E13<S<2X< 1E13.S< 

S2)I<IE13.5<2X<1E13<S</1 
GO  TO  SO 

4S  Cfl<CMPLX{0<0<0<0 1 

•RITE! 3<101 1 I<CC<CO<CO 
SO  CONTINUE 

C  CALL  MAGCHA  SUBROUTINE  TO  CALCULATE  THE  MAGNETIC  CHARGE  DISTRIBUTION 
CALL  MA6CHA1CV<CATN00<NCCNN<NBCUN0<NNC0ES<NE0GES< 

SNFACES<NUNKNS< ITRAK) 

C  COMPUTE  THE  TRANSMISSION  AREA  ANO  THE  TRANSMISSION  COEFFICIENT 

CTS<CMPLX<C<0<C.O) 

00  100  I<1<NUNKNS 
100  CTS«CTS*CV(I1«CCNJG<CI<I}) 

TS«REAL<CTSI/<2<0*AIMP| 

ARE AINbAREAT*COS< THETA )*ALAM0A*«2 

T-«TS/AREAIN 

«RITEt3<llC)  TS<T 

110  FORMAT<*1*,/////<10X<*TRANSMISSION  AREA  ■• < 1E13<S<2X<* SO#  METERS*< 
S///<I0X<*TRANSM1SSI0N  COEFFICIENT  <*<1EI3<S) 
TCHA«TS/<AREAT*ALAM0A**2) 

■RITE<3<120)  TCHA 
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FOFMAT  (//«10X ••TCHA  «**1E13.S) 

IF{  IRCSaEC .0 )  GC  TO  AG9 

C  SINCE  IRCS  IS  NOT  ECUAL  TO  ZERO.COHPUTE  THE  TRANSMISSION 
C  CRCSS  SECTION 

IFC IRCS.EQ. 1 >  GO  TO  299 
IF! IRCS«EQ«2 I  GO  TO  399 
GC  TO  499 

299  HMTsCMPLXI l.O.O.O) 

HMPsCMPLXI 0.0, 0.01 
GC  TO  450 

399  HMTsCMPLXIC.O.O.O) 

HMP«CMPLX| 1.0,0.01 

450  CAtL  TRANCSfOATNCC«NCONN,NeCUNO.NNOOES.NECGES,NFACES. 
SNUNKNSvCV, 1 TRAK.CIMl 
499  CONTINUE 
RETURN 
ENO 
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SUBRQUT IKE  IKCATAI OATKOOvKCCNNvNNCOESaNEOCES ) 

C  THIS  SUBROUTINE  REACS  TtiC  SETS  OF  INPUT  CATA  ANO  ARRANCCS  THEM 
C  IN  NUMERICAL  OROER.THE  FIRST  SET  OF  DATA  CONTAINS  NQOE  NUMBERS 
C  ANO  THEIR  CCCRCINATESaEACH  NOCE  ALONG  WITH  IT*S  THREE  COORDINATES 
C  IS  READ  ANO  STORED  IN  THE  MATRIX  OATNOO. 

C  THE  SECCNO  SET  OF  DATA  CONTAINS  EDGE  NUMBERS  ANO 

C  THE  NCOES  TO  MHICH  THIS  PARTICULAR  EDGE  IS  CONNECTED.  THIS 
C  INFORMATICN  IS  STCREC  IN  THE  MATRIX  NCONN. 

DIMENSION  DATNOOINNOOES.S) 

INTEGER  NCCNNiNECGES.a) 

CO  to  Isl.NNCCES 
REAOfl.SI  NODE.X.V 
S  FCRMATf la.SFB.S) 

ANsFLCAT  INODE) 

OATNCOINCOE* 1 )=AN 
CATNCOfNCCE.2  JsX 
CATNCDINCOE.3  )=Y 
10  CONTINUE 

DO  20  I-l.NECGES 
REAOIt.IS)  NE.NF.NT 
IS  FQRMATI313I 

NCONNINE.l)=NE 

NCCNN<NE.2)sNF 

NC0NN(NE.3)3NT 

20  CONTINUE 
MRiTEca.tai 

IE  FCRMATI*!*) 

MRITEO.  19) 

19  FQRMAT(20X.*VERTEX  CCCROlNATE  LIST*) 

kRlTE<3.21 ) 

21  FORMAT</20X« •all  DIMENSIONS  ARE  IN  METERS*) 
hRlTE(3«22) 

22  F0RMAT(//1X.*VERTEX  NUMBER* .SX.* X-CCCROINATE* .SX. 
••Y-COOROINATE* I 

DC  30  I«1.NNCCES 
lOUMMYslFIXDATNCOI  1.1  )  ) 

liRITE(3.23)IOUMMV.OATNaO<  l•2)•0ATN0DCI•^) 

23  FORMAT I/3X. I3.11X.1E13.5.2X.1E13.S) 

30  continue 
MRiTEo.ae) 

28  FCRMATI*!*) 

IIRITE<3«2S) 

29  FORMATI/10X.*EOGC>VERTEX  CONNECTION  LIST*) 

DC  40  I«1.NE0GES 

WRITE (3.31 INCCNN I I.l >•^CC^N(I•2)•NCCN^f 1.3) 

31  F0RMATI/3X.*E0GE*.I3.1X.* IS  CONNECTED  FROM  VERTEX* vlX. 13. IX. 
S*TC  VERTEX*.1X. 13) 

40  CONTINUE 

RETURN 
END 


Ji 
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SUBROUTINE  6EaM< NCONN • NBCUNO • I TR AK • 1 M IN • NEOGES*NNOOES« 
SKSE«NAPTRE.N»-ANCL  I 

C  THIS  SUBROUTINE  USES  THE  INPUT  CATA  TO  POPP  TRIANGULAR  PATCHES* 

C  THE  INFORMATION  CONCERNING  THE  TRIANGLE  AND  ITS  ASSOCIATED 

C  EDGES  IS  STORED  IN  THE  MATRIX  N80UN0.  ITRAK  AND  IMIN  ARC  T«0 

C  AUXILIARY  VECTORS  NEEDED  IN  THE  PROGRAM* IN  CASE  OF  AN  OPEN  BOOT* 
C  IT  CALCULATES  NUMBER  OF  APERTURES  AND  NUMBER  OF  HANDLES* 

C  IT  ALSO  LISTS  THE  SURFACE  EDGES  ASSOCIATED  MITH  EACH  APERTURE* 

INTEGER  NC0NNINE0GES*3I*NGCUN0(1S0*4)*ITRAK(NE0GES) 

INTEGER  IMININEOGESI 

CCMMCN/IF/IFACe 

IFACE^O 

NFtsO 

NFZsO 

DC  too  IJst.NECGES 
ICQUNTaO 
NlsNCONNI  IJ*2  ) 

NZaNCCNNI 1J*3  > 

CO  10  I«1*NE0GES 
CO  10  J=2*3 
IFIl.EC.lJI  CC  TO  to 
NAaNCONNII*J) 

lF(NA.Ea*Nl*aR*NA*Ea*N2)  GO  TO  6 
CC  TC  10 

6  ICOUNTsICOUNTAl 

ITRAKI lCaUNT)=l 
10  CONTINUE 

MARK1«0 
MARK2S0 

75  CONTINUE 

Xl»I 

llalTRAKIKl) 

CC  IS  |s2*ICCUNT 
IFIITRAKf II.LT*It I  GO  TO  12 
GO  TC  IS 

12  llsITRAKfl) 

XlsI 

IS  CONTINUE 

IFIMARKI .EC* ICOUNTI  GO  TO  tOO 
IFIIl.GT.IJ)  CC  TO  20 
GO  TC  31 

20  CCNTINUE 
N3«NC0NN< 11*2} 

NAsNCONNC 11*3  I 

IF(N3*Ea*M  .CR.NS.EC.NB)  GO  TO  21 
IF(N4*Ea*M*CR.N«*E0*N2}  CO  TO  22 

21  NE*>N4 

GC  TC  23 

22  NB«N3 

23  CONTINUE 
ICC«0 

CC  2S  I«1*NEC€ES 
CO  2S  J«2*3 
IFIl.EC.Il)  GC  TO  25 
NCsNCONNIItJ} 

IFINC.EO.NB)  GO  TO  24 
GC  TC  25 

24  ICOalCCAl 
IMINIICO)*! 

CCNTINUE 
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DC  30  1^1, ICO 
lAsIMlMl) 

IF(NI•E<3.KCC^^(IA•2l•Ofl>^t•EQ«NCONNfIA•3))  €C  TO  29 
lF<^^•EO•^CO^^<IA•^l•OR.^^•Ea•NCCNN(fA•^||  GO  TO  29 
€C  TC  30 

29  l2sIA 

GO  TC  32 

30  continue 

31  continue 

ITRAKI  K1  )«NEOGE£4-l 
MARK1sMARK1«1 
GO  TC  79 

32  IF<12.Lt.lJ>  GO  TC  74 

IFIIFACE.EO.OI  GO  TO  39 
NFlsNBCUN0<IFACE.21 
NF2=NeOUNO(IFACE«3l 

3S  IFI lJ.Ea.NFt.AN0.I2«EQ.NF2l  GO  TC  74 

lFACEaIFACE«l 
NBOUNOI IF ACE • 1 )= IF ACE 
NBOUNOIlFACEaSMIJ 
NaCUN0<IFACE.3l=Il 
NBOUNOI IFACE .4 }s 12 
|IARK2sMARK2»  1 
MARK|3BMARKt#  1 
IF(RARK2.EC.21  CO  TC  100 
ITRAKIKl I-NE0GES41 
CO  TO  75 
74  CONTINUE 

ITRAKIKl }»NE0GES4t 
MARKt«MARKt«t 
C-C  TC  75 
100  CONTINUE 

NSEsO 

OC  120  IsUIFACE 
00  120  JS2.4 
ISEOGEaNBCUNOII.Jl 
NCCUNTaO 

OC  125  Ksl.lFACE 
00  125  M«2«4 

IFI J.EO*K«ANO.J.EO,M|  GO  TO  129 
IF 1 1 SEOGE«EQ«NBCUNO 1 K. F 1 1  NCOUNTeNCOUNT4 1 
125  CONTINUE 

IFiNCOUNT.EQaOl  GO  TO  119 
GO  TO  120 

119  NSEsNSEAl 

ITR AK INSE Is I SEDGE 

120  CONTINUE 
NAPTRE»0 

IFCNSE«E0«0I  go  to  991 

00  147  K|sl«NSE 

IlsITRAKf K1 I 

00  145  3sKl*NSE 

IFIJ.EO.KII  CC  TO  145 

IFIITRAKI JI.GT.lTRAKCNIIl  GC  TO  145 

lOUMMVs ITRAKIKl  I 

1TRAKIRI1«ITRAKIJ) 

ITRAKI4)*10LFNY 
149  CONTINUE 

147  CONTINUE 

00  199  l«l.NSe 


15« 


ININ(  llslTfillKdl 
CCHTIKUE 
df^lTEf  3*  15<«> 

1599  FCBM/IT<*  1"*//*I0X,  •BOUNOARY  CONTCUR  LIST'I 
t6t  CC  900  Jx|«MS£ 

IFIIMlNf J>«AE«0)  GC  TC  109 
.';0  TO  900 
169  IJxIRINfJl 

IRlNf  J}xO 
NAPTRExNAPTREAl 
liRtTE(3«t601  I  NAPTRE 

1601  F0RI»AT<//1X,  •APERTUPe».I3«2X**C0ASISTS  OF  THE 
S  FCLLOaiAG  BCUACARY  ECGES'I 

»RITEf3«16021  IJ 

1602  F0RRATtlSX«I3) 

NtxNCOKNl 1J.2 1 
R2*NC0hN( 1J«3 I 

ITS  DC  600  Kxl•^SE 

fF(IMIN<K>«NE«0)  GO  TO  1T9 
eC  TC  600 
179  IKxlPlNfKl 

RK  |xNCaNN<lK«2) 

FK2sRCCRR<fK«31 
IF! N2«Ea-hKl  )  GC  TC  190 
IFCN2.Ea«NK2 1  GO  TO  19S 
60  TC  600 

190  ||RIT£I3«  1603)  IK 

1603  FCRPAT<1SX«13I 
IMINlKlsC 
N2«NK2 

|F<R2.EC.MI  GC  TO  900 
GC  TO  600 

195  «R1TE<3« 16031  IK 

IKINiKlxO 
R2XRK1 

IF(R2«EC«M)  GO  TO  900 
600  CCRTIRUE 

1F(N2«R6.MI  GC  TO  ITS 

900  CCRTIRUE 

00  901  |X|«RSE 

IFf IKINC 1 1*RE«0I  GO  TC  161 

901  CONTINUE 

991  NHANOLxl>< IFACE-NEOGES^NNCOESANAPTRE 1/2 

RETURN 
END 
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SUeROUT INE  AJU^C(NJUNC•^BaU^O•NEOCES•NFACES) 

C  THIS  SUBROUTIKE  TAEULIZES  ThE  JUNCTICNS. 

INTEGER  NJUNC<NE0CES.3}.N80CN0< 150*4 ) 

CC  10  l*l*hEOGES 
NJUNCll (llxl 
NJUNCt l*2|x0 
NJUNC(1.3)xO 

10  CONTINUE 

00  70  Ixl.NFACES 
NS^NeOUNOl 1.2) 

N3sNe0UN0< 1*31 
N4sN80UN0f 1*41 

IFCNJUNC(N2.2I.EQ*0I  GO  TC  20 
N JUNC<N2*3 l3| 

GO  TC  30 

20  NJUNCCN2*2  }s| 

30  IF<NJUNC(N3*2)*E0.01  GC  TO  40 

N JCNC(N3*3}xI 
GO  TO  SO 

40  N JUNCf N3*2 }=1 

SO  1F|NJUNC{N4*2)*EC.0)  GC  TC  «0 

NJUNC(N4.3M1 
GC  TC  70 

60  NJUNC(N4.21=I 

70  CONTINUE 

hRlTE(3.80> 

80  FCRNATl*!* •20X**JUNCT10N  LIST**//) 

00  too  |3l,NE0GES 
||RITE<3*90)  INJUNCl  l.J)*J»t*3> 

SO  F0RliAT(/3X**E0GE*.13.1X.«  IS  THE  JUNCTION  OF  FACE  *  •  1 X.  13*  IX*  •AND  FA 

SCE* (tXvIS) 

100  CONTINUE 

RETURN 
ENO 


SUBROUTINE  CUROIRf NCONN«N8CUNO.NEACES*NE06ES.1M|N«NSE) 

C  IN  THIS  SUROUTIKE  .NORMAL  VECTOR  TO  THE  SURFACE  IS  CALCULATED* 

C  THE  NCRMAL  VECTOR  IS  CfiTAINEC  EY  LISTING  THE  ECGES  ASSOCIATED 
C  VITH  EACH  TRIANGLE  IN  A  SEQUENTIAL  MANNER.HERE  THE  USER  HAS 
C  THE  CHCICE  OP  SELECTING  THE  DIRECTION  OF  NCRMAL.BUT  IN  THIS  CASE 
C  THE  USER  SHOULC  SUPPLY  THE  INFCRMATION  IN  A  PRESCRIBED  MANNER* 

C  FOR  DETAILS.PLEASE  REFER  TO  THE  REFERENCE  SITED  IN  THE  NOTE* 
INTEGER  NCCNN(NEOGCS.3I.NBOUNC( ISO* 4  I. IPINI NEOGES > 

INTEGER  IMAXI6) 

Iltal 

IMINIIMIat 

MlsO 

DO  S9S  iJKsl.NFACES 
IKS  UK 
OC  2  Is 1*6 

IMAXIIIsQ 

2  CONTINUE 

IFLAGsO 
OC  4  J«1*IM 

IFdJK.EQ.lMINC  Jll  GO  TO  I 
IFLAGsI 
A  CONTINUE 

IFIIFLAG.EQ.il  60  TO  9GS 

I  JIsO 
llaO 
LlaO 

I2sNeOUNO< 1K.2 ) 

I3aNeCUN0(  IK.3  I 
I4aNBOUNO( IK. A) 

NlsNCONN<I2*2) 

N2sNC0NN<I2.3l 
NSsNCGNNf 13.21 

IF(N3*Ea*Nl.CR*N3*EQ*N2>  GO  TO  S 
GC  TC  10 

5  N3sNCONNll3*3) 
to  CONTINUE 
11*1 

ITENP*I2 

II  DO  20  IJst.NFACES 
DO  12  Jsl.lM 

IF(1J*E0«IMIN(J1)  GO  TO  20 
12  CONTINUE 

J2>NB0UN0(  U.2  1 
33aNe0UN0llJ*31 
JA*NBOUNO< 1J.4I 

1FC1TEMP*EQ*32.0R*ITEMP«EO*33*OR*ITEMP*EO.J4|  60  TC  IS 
GO  TC  20 
1 S  IL* 1 J 

GO  TO  2S 

20  CCNTINUE 

IFlIl.EQ.I.ANO.Jl.EO.l I  GO  TC  21 
Jlsl 

ITEPP«I3 
60  TO  II 

21  |FIIl*EQ*l*ANO*Jl*EO*l*ANO*Lt*EQ.ll  GO  TO  23 
Ll«l 

ITEMPalA 
GO  TO  11 

23  lFfMI*EO*l)  GC  TO  009 
MIsRIAl 
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iKi:  UK 
GC  TO  1 

25  KNlsKCCNN< 1TEMP.2) 

KN2sNC0NN< ITEMP.5I 

IFCni .EO.KNl.OK.M .EG.Kh2)  CO  TO  35 

KN3SN1 

60  TO  40 

35  IFCK2«EC.KN1.0P.K2.EC.KN2>  CO  TO  35 
Kh3ah2 

GO  TO  40 

36  KK3sK3 

40  J2=KBCUK0( 1L*2 ) 

J3sN80UN0(  1I..3} 

J4«KBCUK0(IL«4) 

IF(J2*E(3«ITEMPI  CC  TO  59 

IFCNCONNt J2«2).EQ.KN1.0R«NCCNM J2»2>.E0.KK2»  60  TO  57 
KK4SKCCFKC J2.2I 
60  TC  62 

57  KN43:NC0NN<  J2«3) 

6C  TC  66 

59  IFIKC0NNfJ3«2)  .EO.KM  .CRaKCCKM  J3.2  >.Ea.KK2>  GC  TO  61 
KK4sKC0NN( J3«2> 

GC  TC  66 

61  KK4=KC0NN< J3.3I 

66  COKTINOE 

IFllM.EO.l)  GO  TO  115 
IF(1TEMP«EC«IN^X<6I}  CC  TO  109 
IMAXl l)sIM4X<S) 
lMAX(2lsIMAXf 4  I 
IMAX<3>slMAX(6} 

GC  TC  Its 

109  IMAXC  t)s|>lAX<4  ) 
lMAX(21sINAX(6} 

IWAXOlsIKAXiSI 
US  IF(Mt.NE«l)  GO  TO  175 

1F<ITEMP«EC«KBCUKC<IJK«4)>  GO  TO  165 
I MAXI  1 lahBCUKCf IJK.4) 

1MAX(2)=NBQUN0<I JK.21 
1MAX(3)  =  NBCUKC(UK.3) 

MlsO 

GO  TO  175 

165  IMAXl 1 )=KBCUKO( IJK.31 
lMAX(2)=KeCUK0<IJK»41 
lMAX€3>aKBCUK0(lJK.2) 

Ml=0 

1 75  KOUMMVsKNJ 

00  too  1>1<2 

IFCI.EQ.l.AKO. IM.NE.l 1  GO  TO  99 
103l«( I-l 142 

IFI 1TEMP.EQ.12}  GC  TO  79 

IF! ITEMP.E0.13 I  GO  TO  89 

IF(M  •Ea«KN3.AKC.K2.EG.KKl  1  GC  TO  69 

fF(Ml«EQ.KM.AK0.K2.EG«KK3l  CC  TC  69 

IMAXf ICI«13 

IMAXI I042)sI2 

CO  TO  99 

69  IMAXf  I0>sl2 

IMAX<IC«2l3l3 
GO  TO  99 

79  MNl«NCONNf I3«2f 


^^2=^C0^^< 13,3 ) 

lF(NM.EQ.KM.AN0«NN2.EaaKK3>  60  TC  El 
IF(NKt.EC*KN3«ANC«NN2«EC,KM  I  60  TO  61 
IMA)((I0I=14 
1MAX( I0»2)=I3 
6C  TC  99 

61  1FAX(IC)=I3 

IMAXf  tO«2l=I4 
60  TO  99 

89  IF(Nl.EC.Kh3.Ar«0.F2.EG>KM)  CC  TC  91 
IF(Nt.EQ.KNl.ANO.K2«EQ.KN3)  60  TO  91 
IFAXl 101=14 
IMAXC 1042)=I2 
GO  TO  99 

91  IMAXCI0>=I2 

IMAXl 10421=14 

99  KN3=KN4 
12=  J2 
13=  J3 
I4=J4 

100  CCNTINUE 
Kh3=K0UFMV 

NA1=NC0NN(  IMAXC 11,21 
^A2=^CC^^ ( IMAX C 1 ),3} 

^el=^C0^^( 1MAX<4 ) ,2) 

NB2=NC0NN{ 1MAX(41.31 

IF(Mei.Ea«KAl.CF.Nei.E0.NA21  GO  TO  12S 
lFIMa2.EQ.NAl.CF,Ne2.EQ.KA2 1  GC  TC  12S 
101IMMY=1MAX(6  1 
1MAX(6 )=IMAX(4 1 
|MAXC41=I0UMMY 
12S  IMAX(21sITEMP 
1MAX<S)=ITEMP 
IFCIM.ME.ll  CO  TC  149 
NBOUNOC IK,21=1MAX(11 
NECUNOC  IK, 31= IMAXC 2  I 
KeCUN0(IK,4}=IMAX(31 
149  MSCONOt  IL,21=1MAX<61 
NEOUNOC  1L,3I=1MAXCS) 

NBCUMOC  IL ,41  =  11) AX <41 

1M=IM41 

IMINC  IMMIL 

1K=IL 

IFC IM.EO.NFACESl  CO  TC  1000 
GO  TO  1 

999  CCNTIKUE 

1000  CONTINUE 
IFCNSE.EO.Ol  GO  TO  1001 
kRlTEC3,961 

96  FOFMATC*!*! 

MRITEC3,1021 

102  FCRMAT<10X.*LIST  OF  EOGES  ANO  VERTICES  BOUNDING  EACH  FACE*1 
00  1999  1JX=1,NFACES 
12=NeOUNO( IJK.2  1 
ISsNBCUNOC IJK,3  1 
14=NB0UN0< IJK*4  1 

IFCNCONNl 12,21,Ea.NCONN( 13,211  GO  TO  1005 
IF<NCCNNC12, 21. EO.NCONNC  13,311  GO  TC  1005 
N1=NCCNN< 12,31 
60  TO  1006 
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1005 

1006 


1010 

1011 


10  IS 
1016 

1050 

1099 

1001 


Ms^CC^^(  12*21 

IFINCONNI  1^*21. EO.NCONM  14*21  1  GO  TO 
IFf ^C0Nh( 13*2  l.EC.NCONNl 14*3  I  1  60  TO 
^^-^CO^^( 13 .3 1 
GO  TO  toil 
N2=^C0^^( 13*21 

IFCNCONM  14*2  1  .EO.NCONM  12*2  1 1  GO  TO 

1F(NCONN(14*21«EO.NC01«N(12*3I1  GC  TO 

N3=NCONK< 14*31 

GO  TO  1016 

N3SNC0KM  14*21 

CC^T1^UE 

WR1TE(3*10S01  UK* 12* 13* 14*Rt«N2*K3 

F0RMAT1/3X**FACE*  *13.1X** IS  BCUNOEO 

IX* 13*2X* • ANO  VERTICES** lX«13*tX*13* 

CCRTINUE 

RETURN 

ENC 


1010 

1010 


1015 

lOlS 


EY  E0GES**1X*I3*1X 
1X*131 


* 


13* 


110 

111 


107 


s 

6 


10 
1  1 


15 

16 


SUEFQUTINE  BCOPARCOATNCO.NCCNK.NeOUNO.NhCOCSaNEOGESaNFACES* 
SNUNKNS«NSE»NAPTAE«N>-ANOt.  ) 

01  PEAS  I  Ch  0ATNCC(ANQ0ES«ai.AL€3l«»><3>«ReH(3) 

INTEGER  NCCNN<NE0CES*3),KeCi;K0(  ISO .4).  IS  131 
COMMON/ VOL / VOL  LME 

CCPMCN/MOOl/AREAT.RAV* IKMAX* IKM IN .RMAX. RMI N 
CCMM0N/M002/AVAREA* JMAX.ARMAX.JMIN.ARMlAvKMlNsRATIO 
HRITECJ*  1101 

FCRMATI •1*>///2SX.*800V  PARAMETER  LtST*«///l 

•RITE I3.1 1 1 1  NNCOES.NECCESaNFACES.KSE.NUNKNS.NAPTREvNHANOI. 

FORMATC/lOXa'NOMOER  OF  VERT ICE S** • 1 X • 13 • 

S//10X««NUMeER  OF  E0GESs*«lX«13, 

CF  FACES** •IX«t3« 

OF  eOUNOARY  E0GES3*«1X«I3* 

CF  INTERIOR  ECGESCNO.OF  UNKNOHNSI  **»1X«I3« 

CF  APERTURESIEOUNDARV  CCNTCURSI  «*«1X«I3« 

OF  HANDLES  **«1X.I31 


S//tOX«*NUMeER 
S//IOX* •NUMBER 
S//10X.* NUMBER 
S//IOX«*NUMeER 
S//10X. •NUMBER 
AREAT*0.0 


ALTsQ.O 

00  199  IJK*1»NFACES 

ISsNBCiiNOI  UK. 2  I 

I3*NBCUND<I JK.3I 

14*N60UN0< 1 JK.4 } 

lS<t  1*14 

1S<2>*12 

1501*13 

IFINCONNI 12.2 I.EG.NCONNl I3.2I1 
IFlNCONNl 12.2).EQ .NCONNl  13.3  I  1 
N1*NC0NN< 12.3) 

GO  TO  6 

NlaNCCNN<12*2) 

IFINCONNI I3.21.EC.NC0NK(I4«2II 
IFINCONNI t3«2) .EO.NCONNl 14.311 
N2*NC0NN< 13.31 
GO  TC  11 
N2*NCONN< 13.21 

IF|NCCNMI4.21.Ea.NCaNN<  12.211 

IFINCQNNI 14.2 l.EO.NCONNl 12.3  1  1 

N33NCCNNC 14.31 

GO  TC  16 

N3*NCCNN<I4.21 

CONTINUE 

XI*0ATNC0<N1.2} 

Vl*CATNCOf N1 .31 
X2*OATNOO(N2.2} 

Y2*CATN00(N2.31 
X3*CATNC0<N3.2 I 
Y3*0ATNC0f N3.31 

Afl3*<X2-X114<Y3-YI l-CX3-Xll4fY 
AflEA*ABS<AR3 1/2.0 
AREAT*AREAT4AREA 


GO  TO  5 
GO  TO  5 


GO  TC  10 
GO  TC  10 


60  TO  IS 
GO  TC  IS 


2-Vll 


AL<31aSQRT<<X2-XI l•♦2♦^ Y2-YI1442I 
AC<1  1*SCRTUX3-X2|4424€Y3-Ya|44  2| 
AL«21*S0RT«CX1-X314*2*<YI-Y314*21 


ALT*ALT4ALI1  14AL<214AH3> 
H(1 la2.04AREA/AL( 11 
H(21*2.0«AREA/AL<21 
H<31**2.04AREA/AL<3  1 
RBHll 1*H(11/ALC11 
R6H<2|sH<21/ALf21 


PBHCa  »  =  BC3)/AL<2) 

IFCIJK.EQ.ll  GC  tC  19 
00  17  1=1.3 

IF(AL(  I  l.LE.RMIM  GO  TO  19 
17  CC^TI^UE 
GO  TO  21 

19  I•*I^=IS(1) 

ALKIN=A1.<1  1 
OC  20  1=1*2 
AR=AL(  I*ll 

IF(ACFIK*1.E*AR  1  GC  TO  20 
ALFIK=ALI l«l > 

IMIN=ISf 1«1> 

20  CCATir>UE 
FMINsALFIN 

21  IFdJK.CQ.ll  GO  TC  24 
OC  22  1=1*3 

lF<AL(t l.GE.PWAXl  GC  TC  24 

22  COATINUE 
GC  TC  26 

24  IKAXslSlll 
ALMAX=AL< 1 1 
OC  2S  1=1*2 
AP=AtCI41 1 

IF( ALMAX.GE.ARI  GC  TO  25 
A(.AAX=AL(  1411 
1MAX=IS1I41 1 

25  CONTIKOE 
RMAXsALXAX 

26  IF(1JK*EC-11  GC  TC  35 
00  28  1=1*3 

IFlRBhl  D.LE.RATIOI  GO  TO  35 
28  CCATINUE 
GO  TO  4  1 
35  KM1N=IJK 

RATIC1=FBHI1 1 
00  40  1=1*2 

IFlBATI01«LE*RBFf l+lll  GC  TC  40 
RATIOlsOBHl 141 1 

40  CONTINUE 
RATIC=RATIC1 

41  IFIIJX.EO.I)  GC  TC  198 
IFIAREA.LE.ARMINI  GO  TO  191 
GC  TC  192 

191  AHMIN=AREA 
JN1N=1JK 

192  IF<AREA.GE*ARNAXI  GO  TO  196 
GC  TC  199 

196  ARNAX=AREA 
JMAX=  UK 
GC  TC  199 

198  ARXAX«AREA 
JMAX=  UK 
JNIN=IJK 
ARI>IN=AREA 

199  CONTINUE 

RA¥*ALT/FL0AT<NEC6ES4NBCGES-NSe  I 

AVAREA*AREAT/FLCATlNFACeS> 

FACES»»»l*0/AVAReA 

ikxax=ixax 
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nfnin-mth 
1lRIT£(3.20S> 

205  FOWM/lT<//////////2SX,*»OOeHh6  LIST»///I 

llRirE(3.20O  AREAT 

206  F0RMAT(/10X. ■SURFACE  AREA  CF  THE  SC ATTERER* • • 16 1 3«S* I  X. 
S'SO.FETERS* ) 

206  «RITE(3«20S>  RAV • I FAX .RRAX • IN IF .RRIF 

209  FORMAT(/10X,*AVERA6E  E06E  LENGTHS*. IE 13.S • 1X« • METERS* 

S//10X. •MAXIMUM  EDGE  LENCTHIEOCE  NO. • • 13. IX* ■ |s •• IE  13. 

•  •METERS*  • 

S//10X* •MINIMUM  EDGE  LENGTHCEOGE  NO* • *13* 1 X* *  I** • 1 E 13.5* 1 X* 

•  •METERS*  > 

MRITE(3*210)  AVAREA. JMAX.ARMAX.JMIF.ARMIK 

210  FORMAT</tOX**AVERAGE  FACE  AREA  >•• IE I3.S* IX** S0*METERS* • 

•  //tOX«*MAXIMUM  FACE  AREA  IFACE  NO* •• 13* I X • • I* • • IE  13 *5* I X* 

•  •SO.METERS*  * 

•  //tOX* •MINIMUM  FACE  AREA  (FACE  NO* • • 1 3* I X* • )«• • IE  13*5* 1 X* 

•  •SC.METERS*  1 

«|R|TE(3*2I1)  KMIN.RATIC.FACESM 

211  F0RMAT(/10X*<MtN(MUM  FACE  HEIGHT  TO  BASE  RATIO  I 

•FACE  NC***I3*1X*«|s«*|E13*8*//10X**AVERA6E  NUMBER  OF  FACES  PER 

•  SOUARE  METERS* .lElSaSI 
RETURN 

EFC 


CJI  « 


69 


SUBROUTINE  YMATRXICV.OATNCC.NCONN  .KBOUNCcKIkCOESaNEOCeS* 
ShFACES.NUNKNS* ITRAK.Cl  I 

C  IN  This  SUBROUT  INE.ME  COMPUTE  THE  MATRIX  ELEMENTS  AS  DESCRIBED 
C  IN  THE  REFERENCE  SUED  IN  THE  NCTE.FIRST  ME  CALCULATE  THE 
C  POTENTIAL  QUANTITIES  OVER  A  SOURCE  TRIANGULAR  REGICN  AT 
C  THE  CENTRCIO  OF  A  FIELD  TR I ANGLE.THEN  HE  PUT  THIS 
C  QUANTITIES  HITH  APPRCPRIATE  MULTVPLING  FACTCRS  IN  DIFFERENT 
C  ROMS  AND  COLUMNS  OF  V-MATRIX. 

IMPLICIT  CCMPLEX  (C) 

REAL  COS.CABS 

COMPLEX  CVINUNKNSaNUNKNSI •CIINUNKNSI.HTHETA.HPHI 
CCMPLEX  CVTEMPfSO.SOl.CITEMPCSOI 
COMPLEX  CS<3l«ETHETA«EPHl«HX«Hy«H2.H0GTT 
DIMENSION  OATNOO<NNOOES»3}«TMAT(3*2) 

INTEGER  NCCNNINECCES.Sl.NBCUNCf ISO. 4 )• I TRAK CNE06ES I 

COMMCN/KKK/AK.PI 

COMMON/PARAM/ThETA.PHI.lF lELO 

CCMMCN/FIELC/ETHETA.EPhl.ALAMOA 

CCMMCN/OIELEC/CIA.OIB 

NFIEL0=1F1EL0 

ClsCMPLXIl. 0,0.0) 

C2«CMPLX(2.0,a.0) 

2  DO  1040  lOlE-1,2 

OlEsSCRTCOIA) 

1F(I0IE.EQ.2)  OIE=SCRT(Oie) 

C  CALCULATE  ELECTRICAL  PARAMETERS 
PIa3,141S926S 

AKs2.0*PI/ALAM0A/SQRT(0IA>4CI£ 

VGL33,OE40E/01E 

AOMEGAsAK*VEL 

CCNSTlsCMPLXl l«0E-O7,O.0)*CMPLXf0.0,AOMEGA) 
CONST23CMPLX(S.OE409. 0,0 )4CMPLX 10.0, 1«0/ACMEGA)/01E««2 
AlMP=120.0«PI/OtE 
CIMP«CMPLX|AIMP,0.0) 

CONST lsCON<TI/(ClMP*C IMP) 

C0NST2=C0NST2/<CIMP4CIMP) 

C  CALCULATE  HTHETA  ANC  hPHI 
HPHI*ETHETA/CIMP 
HTHETA=-EPh I/C  IMP 
CTl«CMPLX(CCS(ThETA),0,0) 

CT2»CMPLX<SIMTHETA).0«0) 

CPHIlaCMPLXICOSIPHI ),Q.O) 

CPHI2»CMPLX(SlN<PHlf ,0.01 

C  COMPUTE  THE  AMFLITtOE  OF  INCIDENT  MAGNETIC  FIELD  IN 
C  TERMS  OF  HX,HV  ANO  hZ, 

HX3HTHETA4CTl*CFhlt»HPHI4CPH12 
HT<*HTHETA«CT  14CPHI24HPHI4CPH.il 
HZa-HTHETA4Cf2 
NSE«NEOGES'>NUNKNS 
00  999  IJKsl.NFACES 

lFClJK.NE,l.AN0,NFIEL0.6T,t>  GO  TO  999 
C  CETAIN  THE  EDGE  NUMBERS  OF  THE  SOURCE  TRIANGLE, 

12«NBGUN0<IJK,2) 

I3«NB0UN0< IJK«3) 

14«NBCUNC1 IJK,4) 

C  OBTAIN  THE  VERTICES  CF  THE  SCURCE  TRIANGLE, 

IFINCONNI I2*2),EC,NC0NN(13,a))  GO  TO  S 
IFINCONNI I2,2),EC,NCONN(I3,3ll  GO  TO  S 
NI«NC0NNII2*3) 

GC  TC  6 


5  M»NCCNN(  12.2) 

6  IFCKCOhNt  I3.2).EC.NC0NN(I4«21I  CQ  TO  10 

IFCKCCNM  13,2  I.EQ.^Ca^^(  14.3  I  1  CO  TO  10 

K2>^CQ^K( 13.3) 

GO  TO  11 

to  ^^s^COK^{ 13.21 

It  IFCNCONNl I4.2).EC.^CCN^(12•2))  CO  TC  IS 

IFINCONNI 14.2). EQ.NCONNl 12.3)1  CO  TO  IS 

^3^^CC^^( 14.31 
GO  TC  16 

15  N3-NCONM  14.2  ) 

16  CC^TI^UE 

C  DETAIN  THE  COORDINATES  CF  THE  VERTICES  CF  THE  SCURCE 
C  TRIANGLE. 

XlsCATNCOIM.a  ) 

Yl=DATNCO<M  .3) 

X2=OATNOO(N2.2) 

V2^CATNC0(N2.3 1 
X3sCATNCO(N3.2 I 
V3=OATNOO(N3.3) 

C  calculate  the  area  CF  THE  TRIANGLE  BY  TAKING  THE 
C  MAGNITUDE  OF  THE  VECTOR  CROSS  FRCGUCT  CF  THE  TMC  SIOES. 
AR3=<X2-X1)4<Y3-Y1)-I X3-X 1 )•( Y2-V1 ) 
AREA=AeS(AR3)/2.0 

C  DETAIN  THE  LENGTHS  CF  EACH  SIDE. 

R2MR1M=SQRT((X2-X1)«*2H(Y2-Yt)««2) 
R3I«R2M=SQRT(<X3-X2  )*42«<  Y3-Y2  )«42  ) 
fl|KR3HsSaRTC(Xl>X3)«*2«<Vl'-V3|4«2) 

C  CBTAIN  THE  HEIGHTS  OF  EACH  VERTEX  KITH  RESPECT  TC 
C  THE  CORRESPONDING  OPPOSITE  EDGE. 

CHlsCMPLXI 2.04ARCA/R3MR2M.0.0 ) 
CH2sCMPLX(2.0*AREA/RlMR3M«0.OI 
CH3«CMPLX(2.0*AREA/R2MR1M.O.O) 

C  NOW  CALCULATE  THE  PARAMETERS  OF  THE  FIELD  TRIANGLE. 

00  ASS  IJst.NFACES 

C  OBTAIN  THE  EDGES  OF  THE  FIELD  TRIANGLE. 

32«N8CUN0( IJ.2 } 

J3=NB0UN01 IJ.3) 

J4SNB0UNDI 1 J.4) 

C  OBTAIN  THE  VERTICES  OF  THE  FIELD  TRIANGLE. 

IFINCQNN(J2.2).EQ.NCGNNCJ3«2))  GO  TC  250 

IFINCONNC J2.2).Ea.NCaNN<33.31)  GO  TO  250 

NJlsNC0NNfJ2.3 } 

GO  TC  255 

250  NJ1sNCONN<J2«2) 

255  IFiNCONNf J3.2).E0.NC0NN(J4.2 ) 1  GO  TO  256 

1FCNC0NN<33.2).EC.NC0NN(J4.3))  GO  TC  256 

NJ2sNCCNNf  J3.3) 

GC  TC  2Sa 

256  N J2sNCCNN€ J3.2) 

258  IFiNCONNf J4.2).Ea.NCONN<J2.2))  GO  TC  25S 

IFINCCNNf J4.2).E0.NCaNNf J2.3I)  GO  TO  2SS 

N j3sNCONN< J4.3) 

GO  TO  260 

259  NJ3sNC0NN<J4.2) 

260  CONTINUE 

C  OBTAIN  THE  CENTROID  OF  THE  FIELD  TRIANGLE. 

Xaf OATNOOIN31 .2 ) 40 ATNOOC N J2.2 140 ATNOC IN J3. 2  1 1/3.0 
V«<0ATN00INJ1.3>40ATNC0(NJ2.3140ATNCDINJ3.31 1/3.0 
C  CALCULATE  COMPCNENTS  CF  THE  TESTING  VECTOR 
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C  CCfiRESPONOlNG  TO  EACH  SIDE. 

TMATCl . 1 )= (0ATHCC(KJ2.2 ) 40 ATNCO (N J3 • 2 n/2 • 0-A 
TMAT(1.2>=(CAThC0(KJ2«3 )»CATKC0CNJ3.3) >/2.0>Y 
TMAT(2*  1  )=(OATKCOCN J3«2)*OATKCO|NJt«2> 1/2.0-X 
TI*AT<2«2l=(0ATNCCf  NJ3«3)»CATNOO(NJl«31  )/2.0-Y 
TMAT(3« 1)=<OATACC(KJ1.2}AOATKCO(NJ2.2I1/2«0-X 
TMAT(3•2)s<0AT^C0(^ J1.3)«0ATKC0«NJ2«3I )/2.0-Y 
C  OBTAIN  THE  COORDINATES  OF  THE  FIELD  TRIANGLE. 

X JIsOATNCOIN Jl .2  I 
YJlsOATNCO(NJl«3l 
XW2sCATNC0(NJ2.2 ) 

YJ2s0ATN00fK J2«3) 

X J3sOATNOO(N J3.21 
YJSsOATNCO (NJ3.3 I 

C  DETAIN  THE  LENGTH  CF  THE  EACH  SIDE. 

CS( 1 IsCMPLXISORTf <XJ2-XJ3I««2«<YJ2-VJ3>««2>«0.0) 
CS(2)sCMPLX(SCRT((XJ3-XJI)*«2«(VJ3-VJll««2).0.0l 
CS(3I=CNPLX<SCRT(<X Jl-X J2t««2H<YJl-YJ2>*«2l»0.0) 

IFINF lELO.NE. 1 1  GC  TO  261 
C  CALL  THE  SUBROUTINES  TO  COMPUTE  THE  INTEGRALS. 

CALL  SCAINTIXl •Y1.X2.Y2.X3.Y3.X.Y.CPH1. AREA) 

CALL  VECINTI XI .Yl.Xa.YE.XS.VS.X.Y. 

SCAXSI.CAETA.AREA  ) 

261  IVsO 

C  COMPUTE  THE  VECTOR  ANO  SCALAR  PCTENTIALS  ASSOCIATED  MITH 
C  EACH  EDGE  OF  THE  SOURCE  TRIANGLE. IF  ANY  OF  THE  THREE  EDGES 
C  IS  A  BOUNDARY  ECGE.THEN  THE  CURRENT  COEFFICIENT  OF  THIS  EDGE 
C  IS  ZERO  ANO  HENCE  THE  POINTER  JUMPS  CUT  CF  THE  LOOP. 

C  HERE  CAX.CAY  ANO  CAZ  ARE  THE  VECTOR  POTENTIALS 

C  IN  THE  X.Y  AND  Z-OIRECTIONS  RESPECTIVELY. 

00  4«0  IKsl.3 
IF(IK.EQ.l)  Il»14 
IF(IK.EC.2)  IlsI2 
IF(IK.EG.3)  11=13 
K1=0 

IF(NSE.EC.O)  GC  TO  286 
DO  28S  Jal.NSE 

IFIIl.EO.ITRAKf J))  GO  TO  4«0 
285  CONTINUE 

DO  287  K=1.NSE 

IFC Il.GT.ITRAK(K) )  GO  TC  266 
GC  TC  268 
266  KlsKlNl 

287  CONTINUE 

288  IFCIK.Ea.2)  GC  TO  300 
IFIIK.EO.l)  GC  TC  310 
CFLAG«C1 

IFCNCONNI 11.2 ).EQ.N2.AN0.NC0NN< I l.SI.EQ.Nl)  CFLAG»-CI 
CAXsCFLAG4(CMPLXlXl-X3.a.0)«CPHl«CMPLXf X2-XI.0.0)*CAXSI 
S«CMPLX<X3-X1.0.C)*CAETAI/CH3 
CAYaCFLAG*<CMPLXIYl-Y3.0.0)4CPHI4CMPLX| V2>V 1.0.0 IMCAXSI 
S«CMPLX(Y3-Y1 .O.OIACAETA )/CH3 
CSP0T«CFLAG4CCNST24CPMI«C2/CH3 
CO  TO  37S 
300  CFLA6«C1 

IFCNCONNl I1.2>.E0.N1.AN0.NCCNN(II.3).E0.N3)  CFLAG«-C1 
CAX«CFLAG*(CMPLXCX1-X2.0.0)4CPH|4CMPLXCX2-XI.O.O)4CAXSI 
S4CMPLX I X3-X I .0 .0 ) ACAET A )/CH2 
CAYsCFLAG*<CMPLX(YI-Y2.0.0)4CPH|4CMPLXCY2>V1.0.0)OCAXSt 
S4CMPLX(Y3-YI.0.0I4CAETA)/CH2 
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CSP0TsCFLAG*CaK£T2*CPHl*C2/CN2 
CO  TO  275 
310  CF1.AG=C1 

1F(NC0NN< II.2 1 .Ea«N3.AK0«NC0NM(lt«3)*E0*F2}  CFLACs-Ct 
CAXsCFLAG«<CMPLXOl2-Xl«0.0)*CAXSIACFPEAIX3-Xl«0.0>*CAETA)/CFt 
CAVsCFLAG*(CMPl.X(Y2-Vl«0>0  IPCAXS  I4CMPLXC  V3- V 1  •  0*  0  ) *CAETA | /CH 1 
CSP0T«CFLAG*C0NST2«CPHI*C2/C»41 
375  CONTINUE 
IVslV«l 

C  COMPUTE  THE  TESTING  VECTOR  ASSCCIATEO  MITH  EACH  EDGE. 

C  AGAIN  IF  ANY  OF  THE  EOGES  OF  FIELD  TRIANGLE  IS  A  BOUNDARY 
C  EDGE  THEN  THE  POINTER  JUMPS  OUT  OF  THE  LOOP. 

OC  4S0  IRst.S 
IFdR.EQ.ll  J1=J4 
IF<IP.E0.2I  J13J2 
IFflR.Ea.3)  J1=J3 
LlaO 

IFINSE.EG.Ol  GC  TC  405 
OC  300  Jal.NSE 
IFIJI.EO.ITRAKI Jll  GO  TC  450 
390  CONTINUE 

DO  399  Kal.NSE 

IFIJl.GT.ITRAKIKll  GC  TO  397 
GO  TO  409 
397  LlsLIHl 
399  CONTINL' 

405  crisCMPLx^  MATI  IR.l I.O.OI 
CT2«CMPLX(TMAT  <  IP«2  >.0.0  1 

C  COMPUTE  THE  DOT  PPCOUCT  BETMEEN  THE  VECTOR  PCTENTIAL 
C  AND  THE  TESTING  VECTOR. 

CTEMPsCCNST 1 4 (CAXACT I «CAr ACTE  > 

ARGMNTsX*SIN<THETA)*C0S(PHl)4V*SlNfTHETA)4SlNfPHI) 

C  ARGsCMPL X  (  0. 0  ••>AK  4ARGMNT  ) 

HOCTTsHX ACT 14HV ACTE 
CHTEMPsHOCTTACEXPICAPGI 
IFIIR.EQ.I)  GO  TO  420 
1F(IR. EC. 21  GC  TC  430 
CFLAGsCl 

IFINCONNC J1.21.EC«NJ2.AN0.NC0NN( JI.ai.EC.NJl I  CFLAG»-C1 
GC  TC  440 
420  CFLAG^Cl 

IF(NCONN( J1.2)«E0.NJ3.AN0.NCCNN( J1 .SI.EC.NJE)  CFLAC«-Ct 
GO  TO  440 
430  CFLAGsCl 

IFINCONNI JI.2I.CQ.NJ1 .AN0.NCCNNCJ1«3I«EC.NJ3I  CFLAG«-Ct 
440  IFINFIELO.NE.I  I  GO  TO  442 

CYI JI-LI • I I-K I l»CY <  Jl-L 1 • I I-K I • ♦CFL AG* ( CTEMP-CSPOT  >*CS« IR I 
442  IFIlJK.NE.l.OR.lV.CT.ll  GC  TG  490 

C 1 1 Jl-L 1 J*C I «J l-L 1 1 ♦CFL AG4CSI IR l♦CHTEMP 
490  CONTINUE 
460  CONTINUE 
499  CGNTINUE 
999  CONTINUE 

IFINFlELO.NE.il  GC  TO  1009 
DO  1000  Isl.NUNKNS 
DO  1000  J«1.NUNKNS 
1000  evil.  JloCYl  I.JMCPPLXIO.O.O.OI 

1009  DO  1010  Isl.NUNKNS 

1010  CIIII«C1<1)4CMPLX<2.0«0«0) 

IFIOIA.EO.OIBI  GC  TC  2000 


IF(IOIE.EC«2l  CC  TC  1025 
00  1020  |sl•^U^F^£ 

00  1015  JslaNUNKNS 
CVTEFPl 1 • J l>CY 1 1 • J  I/C2 
1015  CVCl. JlsCOPLXlC.O.O.O) 
CITEMPClIaCKl) 

1020  CIlllsCMPLXIO. 0.0.0) 

60  TO  1040 

1025  00  1035  ls|•NUNK^S 

00  1030  Jsl.NUNKNS 
1030  C«CI.J)aCV<I.J)/C24CVTEI«P<  l.Jl 
1035  Cl(  I)>sCITEMP(l  ) 

1040  COhTINUE 
2000  RETURN 
ENO 


f 

\r 
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SUBROUTINE  VWl  NCC (CY .  O/ITKOCaNCONN. NBOUNO.NNOOES.NEOGES.NEACES.NUNK 
SNSvlTRAK.NJLKC.SICWA.OELO IE •TRICK I 
C  This  SUBROUTINE  CALCULATES  THE  ADDITIONAL  TERM  IN  THE  ADMITTANCE 
C  MATRIX  WHEN  WE  COVER  THE  APERTURE  WITH  A  {LOSSY |  DIELECTRIC  SHEET. 
IMPLICIT  CCMPLEX  (C) 

REAL  COS.CABS 

COMPLEX  CVCNUNKKS.NUNKNSI.CSI3I 
DIMENSION  OATNCO( ANODES .31. TMATC3.2> 

IN TECER  NCONN (NEOCES.3>. NBCUNO flS0.4>.l TRAK ( NEDCE S I 

INTEGER  NJUKC<KE0CES.3I 

COMMON/PARAM/THETA.PHI.IPIELO 

COMMON/FRE/AOMECA 

P1«3.I4IS92«S 

CONST«CMPLXCSIGMA.AQMEeA4CEL0IE/«36«PI*l.0E«09ll 

CONST^CONST •CMPLXI THICK • 0. 0 1 

CIsCMPLXll.O.O.OI 

NFlELOsIPlELO 

IVsO 

NSEsNEOCES-KUKKNS 
DC  999  Isl.NEOCES 
NFlsNJUKCl 1.2) 

KF2«KJUKC(I.3I 
IF(NF2.KE.0I  eC  TC  868 
IV«IV«1 
GO  TO  999 

888  DO  499  M»|.2 

J2aNe0UN0<KFl«2l 
IFIM.E0.2)  J23Ne0UN0{NF2.2l 
JSaKBOUNOIKFl.SI 
1F<M.EQ«2|  J3*KSCLN0(NF2«3I 
J4sNB0UND(NF1.4} 

IF(M.EQ.2I  J4«KeCUN0(KF2.4l 
C  OBTAIN  THE  VERTICES  CF  THE  FIELD  TRIANGLE. 

IFINCONNI J2,2I.Ea.NCQNNlJ3«2ll  GO  TO  2EC 

1FINC0NN(J2.21.EQ.NC0NN( J3.3I }  GO  TO  2S0 

NJ1«NC0NN( J2.3I 
GO  TO  29S 

250  NJ1«NC0NN( J2.2I 

29S  1F(NCONNIJ3.2I.GO«NCQKNC J4.2II  GC  TC  256 

1FCNC0NN(33.2I.EQ.NCQNN(J4.3II  60  TO  25C 

KJ2*NCCNN( J3.3) 

GO  TC  258 

256  NJ2«NCaNN< J3.2I 

258  IFlNCONNf J4.2I.E0.NC0NN(J2.21I  GO  TO  299 

IF(NCONN(J4.2).Ea.NCONNIJ2.3}l  GC  TC  299 

NJ3«NCaNNIJ4«3 I 
GC  TC  260 

299  NJSsNCONNI J4.2) 

260  CONTINUE 

C  OETAIN  THE  CENTROID  OF  THE  FIELD  TRIANGLE. 

XslOATNOOI NJt .SlMDATNCOCNJB.BlMOATNCOIKja.BII/S.O 
V«|  OATNOOl  N  J 1  •  3  >  40ATN00  {  N  J  2  .3 140ATN00  INJ3.3n/3.0 
C  CALCULATE  CDMPCNENTS  CF  THE  TESTING  VECTOR 
C  CORRESPONDING  TC  EACH  SIDE. 

TMATf 1. tl«C0ATNO0(NJ2.2l4OATNCOfNJ3.2ll/8.0-X 
THAT! 1.2 l«COA«NCO(N J2.3 |40ATNC0<NJ3.3) >/2.a>V 
TMATC2.1)«fOATNCO<NJ3.2}4OATNCOiNJ1.2l>/2.0-X 
TMATI2.2>«IOATNOO<NJ3.3l40ATNOOtNJI.3l 1/2.0-V 
TMATC3. I )«(OATNOO(NJ1.2|«OATNOOfNJ2.21>/2.0-X 
TMAT<3.2>«COATNCOf  NJ1.3)40ATNC0<NJ2.3)I/2.0'>V 
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IFfl.E0.J2)  IIV«1 
lFfI.E0.J3)  1IV«2 
IF(I.E0.J4I  llV-3 

C  QETAIM  THE  COCROINAIES  OF  THE  FIELD  TRIANGLE. 
XJls0*TN00tHJl«2) 

VJl«OATNOOfNJl.a) 

XJ2s0ATN0DfNJ2.2) 

VJ2>OATKOOCNJ2.3) 

XJ3sOATNOO(NJ3.2l 

YJ3«CATNC0<NJ3.3) 

C  OBTAIN  THE  LENGTH  CF  THE  EACH  SIDE. 

CSf 1 )«CMPLX(S0RT((XJ2>XJ3)A«2«<VJ2-VJ3)**2>«0.0) 
CSf2)sCMPLX(S0RT((XJ3>XJt )««2F(«J9-VJ1)*«2) .O.O) 
CSI3>sCFPLX(SCRT(<XJl-XJ2)**2F«VJl-VJ2)«*2).0.0) 

AR3sf XJ2-XJt)*tTJ3-VJI)-<XJ3-XJ))*|yj2-VJl I 
ARE A«ABS(AR3  1/2.0 

C  CONFUTE  THE  TESTING  VECTOR  ASSOCIATED  WITH  EACH  EDGE. 

C  AGAIN  IF  ANY  OF  THE  EDGES  OF  FIELD  TRIANGLE  IS  A  BOUNDARY 
C  EDGE  THEN  THE  POINTER  JUMPS  OUT  OF  THE  LOOP. 

DO  450  IRsl.3 
IFClR.EO.l)  JlaJ4 
1F<1R.EQ.2)  J1=J2 
1F<1R.EC.3)  J1*J3 
Ll-O 

IF(NSE.EO.O)  GO  TO  405 
DO  390  Jsi.NSE 

IFl Jt.EC.ITRAKI J) )  GO  TC  450 
390  CONTINUE 

DC  399  Ksf.NSE 
IFIJt.GT.lTRAKfX) )  GO  TC  397 
GO  TO  405 
397  L1SL141 

399  CONTINUE 

405  CTlsCM''iX(TMAT<IR.l  )*TMAT<  IIV.D.O.O) 

CT2sCMPLXfTMAT( IR.2)*TMAT( 1IV.2)«0.0) 

CLCA«CSCIR)*CS< 1IV)/CMPLX<4«AREA«0.0> 

IFClR.EO.l)  GO  TC  420 
IF4IR.E0.2)  GO  TC  430 
CFLAGsCt 

IFCNCONNf  J1.2).EO.NJ2.ANO.NCONNf  JI.3).E0.NJI  )  CFLA6«->CI 
GO  TO  440 
420  CFLAG«C1 

IFfNCONNC J1.2) .EO.NJ3.ANO.NCGNNf J1«3).C0.NJ2)  CFLAG«>C) 
GO  TO  440 
430  CFLAG«C1 

IFCNCONNC JI.2).E0.N Jl.ANO.^CO^NC JI.3I.E0.NJ3)  CFLAG*>C1 
440  IFCNFIELO.NE.l )  60  TO  450 

CYll-IV.Jl-Ll )sCYC l-IV. JI>LI)4CFLAG«C0NST*CL0A*fCTlHCT2) 
450  CONTINUE 
460  CONTINUE 
499  CCNTINUE 
999  CONTINUE 
RETURN 
END 
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SUBROUTINE  MAec»<A(CV.O/>TNOC*NCONN*NaCUNO*NNQOES«NEOGES« 

SNFACESaNUNKKS.lTRAK) 

C  THIS  SUBROUTINE  COMPUTES  THE  MAGNETIC  CHARGE  DISTRIBUTION  ON  THE  APERTURE 
C  AREA«T»>E  CHARGE  DENSITY  IS  COMPUTED  AT  THE  CENTROID  OF 
C  EACH  TRIANGLE. 

IMPLICIT  CCMPLEX  <CI 
REAL  CQS.CABS 
COMPLEX  CVINUNKNSI  .CSOI 
DIMENSION  OATNCOINNODES.3I 

INTEGER  NCCNN(NEGGES.3I«N8CUN0(1S0.A  »• ITRAK (NE06ES > 

COMMON/KKK/AK.Pl 
COMMON/O lELEC/O lA .0  IB 
CISC MPLXI I. 0.0.0) 

VELs3.0E«Oe/SORT<OIB) 

AQMEGA«AK*VEL 

CONSTlsCMPLXfO.O. 1.0/ACMEGA) 

CHARGEaCMPLXIO.O.O.O) 

MRITEia.lOl) 

101  F0RMAT(*1* •///2SX.*SURFACE  MAGNETIC  CHARGE*.///) 
kRITE(3.I02> 

102  FORMAT! IX. 'FACE  NUMBER* • 2SX. •MAGNET IC  CHARGE  DENSITY  ! MEBS/M-M) * ) 
liRITE<3.103) 

103  FORMAT</20X.*REAL*. 11 X.* IMAGINARY* .SX.* MAGN ITUOE* .lOX. *PHASE*  t 

nse»neoges-nunkns 

DO  000  IJMsl.NFACES 
C  OBTAIN  THE  EDGES  OF  THE  TRIANGLE. 
l2>N60UN0f IJK.2 ) 

I3SNBCUN0! I JM.S) 
lAsNBCUNOl I UK .4) 

C  OBTAIN  THE  VERTICES  CONNECTED  TO  THESE  EDGES. 

IFCNCCNN! I2.2).E0.NC0NN! 13.2) )  GO  TO  S 

1FINCONN(I2.2I.EO.NCONNCI3«3II  GO  TC  S 

NlsNCONN! 12.3) 

GC  TC  6 

N|sNC0NN(l2.2l 

IF(NCaNN!13«2).EQ.NCaNN!I4.2))  GO  TO  10 

IFINCONN! 13.2). EC.NCONN! 14.311  GO  TO  10 

N2sNCCNN<l3.3) 

GO  TO  11 

10  N2SNCCNN! 13.2) 

11  IFINCONN! 14.2).EG.NCaNN<I2.2))  GC  TC  IS 

IFfNCDNN! 14.2 I.EO.NCDNN! 12.3))  GO  TO  IS 

N3SNCCNN! 14.3 ) 

GO  TC  16 

15  N3SNCONN! 14.2) 

16  CCNTINUE 

C  COMPUTE  THE  COORDINATES  OF  EACH  VERTEX. 

XlsOATNOOCNI.2) 

VlsCATNOOINI.3) 

X2sOATNOO!N2.2) 

V2>0ATN00!N2«3) 

X3s0ATN00!N3.2) 

V3a0ATN00!N3.3l 

C  CALCULATE  THE  AREA  OF  THE  TRIANGLE. 

AR3sCX2>XI )4!Y3-VI I-! X3-X 1 )«C V2-VI I 
AREA«A8SCAR3)/2.0 

C  CALCULATE  THE  LENGTHS  OF  EACH  SIDE. 

R2MR1MmS0RT!IX2>X1 )4424!Y2>V1 )*«2) 

R3MR2M«SaRT<!X3-X2)*42«(Y3-V2>4*2) 

RIMR3M«S0RTf !X1>X3)**2«!V1«V3)*42)  £ 

\ 

■  ’S 
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CS<1  )sCMPL»(R3*iR2M«0.Q) 

CSt2  )sCI*PLXCRl»«R3R*0.0) 

CSC3)=CMPLX(R2I«R1Ii«0.0) 

C  COMPUTE  THE  MACNET 1C  CHARGE  DENSITY  ON  THE  TRIANGLE. 

CSUMaCMPLX (0.0.0.01 
DO  460  IKal.3 
IFdK.EQ.ll  It«14 
IFIIK.EC.Z)  IlsI2 
IFIIK.EQ.3)  I1«I3 
Kl=0 

IF(NSE.EG.O)  GC  TO  288 
00  285  JsI.NSE 

IFdl.EQ.ITRAKCJll  GO  TO  460 

285  CCNTINUE 

00  287  Ksl.NSE 
1F(I1.GT.(TRAK(K|}  GO  TO  286 
GO  TO  288 

286  K1«K141 

287  CONTINUE 

288  1F(IK.EG.2I  GC  TO  300 
IF(lK.EQ.l)  GO  TO  310 
CFLAG^^Cl 

IFINCONNC  II.21.EQ.N2.AN0.NCCNN<I1.3).E0.N1 1  CFLAC»C1 
GO  TO  375 
300  CFLAG»CI 

(F(NCONN(  11.2). EO.M.ANO.NCCNNt  11.31. E0.N3  I  CFLAGs-Cl 
GO  TO  375 
310  CFLAGsCl 

(F(NC0NN(ll.2).Ea.N3.AN0.NCCNN(Il«3).E0.N2)  CFLAGs-CI 
375  CCNTINUE 

CSUM3CSUM4CFL  AG4CCNST 1 *CS ( IK  I «C V ( 1 1-K 1 > 

460  CONTINUE 

CHCEN>CSUM/CMPL  A ( AREA. 0.0) 

RA1«RCAL(CH0EN) 

RAEsAIMAGICHOEN) 

RAJsCABSfCHOEN) 

RA4>ATAN2<RA2.RA1 ) 

liRITE(3.S01  )  IJK.FAl.RA2.RA3.RA4 

501  FORMAT (2X. |4.eK. 1E13.S.2A, 1E13.S.2X. IEia.S.2X.lEI3.5./) 
CHARGE«CHAFGC4CSUK 

990  CONTINUE 

«RITE(3.502)  CHARGE 

502  FORRAT<////tOX.*TCTAL  MAGNETIC  CHARGE  ON  THE  AREAa  ( •.2E13.S.1X. 
«■  tiEBERS*) 

RETURN 

ENC 
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SUBROUTINE  TRANCS < OATNOO«NCONN.NeoUNO«NNODES«NEDGES*NFACES« 
SNtKKNS»CV»ITRAK«CIM) 

C  IN  THIS  SUBROUTINE.  THE  TRANS  CROSS  SECTION  IS  COMPUTED  AS  A 
C  FUNCTION  OF  SPHERICAL  COORDINATE  ANCLES  THETA  AND  PHI. 

C  PHIl  AND  PHI2  REPRESENT  THE  INITIAL  AND  FINAL  VALUES  OF  PHI 
C  AND  NPHI  REPRESENTS  THE  NUMBER  OF  DIVISIONS  IN  PHI  DIRECTION. 

C  SIMILARLY  THETAl  ANO  THETA2  REPRESENT  THE  INITIAL  AND  FINAL 
C  VALUES  OF  THETA  ANO  NTHETA  REPRESENTS  THE  NUMBER  OF  DIVISIONS  IN 
C  THETA  OIRECTION.THUS.  THIS  SUEROUTINE  COMPUTES  TRANS  CROSS  SECTION 
C  FOR  NPHI  X  NTHETA  VALLES  CF  THETA  ANO  PHI. 

IMPLICIT  COMPLEX  CC) 

REAL  COS.CABS.CCNST 

COMPLEX  CIMINUNffNST.CViNUNKNSI.HMT.HMP 
DIMENSION  OATNOOINNOOES.Sl 

INTEGER  NCCNNINEDGES.a I.NBOUNOC ISO. 4} • ITRAKC NEOGES) 

CCMMCN/FHINCI/FHINC 

COMMON/KKK/AK.PI 

CCMMCN/POLARM/HMT.HMP 

COMMON/OIELEC/ClA.Oia 

READ! I.IGSIPHI l.PH I  2. NPHI .THETAl .THETAa.NTHETA 

195  FORMAT(2F7.2. 13. 2F7. 2.131 
HRITE(3.196I 

196  FORMAT! • 1* .////BOX.* APERTURE  TRANS  CRCSS  SECT lON/SC.H.L. • •///> 
HRITC<3.I97I 

197  FORM AT  CSX. 'THETA (DEGREES )« .SX. •PH I<OEGREES 1 *.SX. 

S3X.*TRANS  CROSS  SECTION  CSC.NETERS/SO.METERS >• } 

0PH13(PH12>PHI1 1/FLOATCNPHIl 

0THETA=!(THETA2-TH£TA1  I/FLOAT! NTHETA  I 

NTHETaNTHETAHI 

NPH^NPHIHI 

OC  999  Ilal.NTHET 

THETA»THETAIHFLCAT<I(-I l♦CTH£TA 

THETAsTHETA*P 1/100.0 

00  998  JJsl.NPH 

PHIsPHI IHFLOATC J J-t I40FHI 

PHI^PHIAPI/teO.C 

ALAMOA»2«PI/AK 

AaMEGAaAK*3.0E«08/S0RT<Ciei 

CONST=< A0MEGA/<364PI) l•l•0E-C«•0I8 

C0NST*CCNST442/<e*PI) 

CALL  MEASURCCIM.THETA.PHl.OATNOO.NCONN.NEOUNO.NNOOES. 
SNEOGES.NFACES.NUNXNS. ITRAK) 

CTCSsCMPLX(0.0.0.07 
OC  499  IJKsl.NUNXNS 
CTCSsCTCSHClMC I JK|4CV(IJKI 
499  CCNTINUE 

TCSsICABSCCTCSI I4424CCNST 
TCS«TCS/ALAMOA4«2 
TCS=TCS/FHfNC442 
THETAOsTHETA*180.0/PI 
PH lO^PH I* 100.0/FI 
«R1TE(3.991)  THETAO.PHIO.TCS 
991  F0RMAT(/3X.1E13.S.4X.1E13.8. 12X.IEI3.S) 

998  CONTINUE 

999  CONTINUE 

1F(CABS(HMT> .EO.FLQAT! Ill  GO  TO  1000 
IFICABSCHMPI.EC.FLOATCIII  GC  TO  1010 
60  TO  1020 
1000  »RITE(3.100SI 

1005  FORMAT4//3X,«TCS/SO.«.L.  OF  THETA  POLARIZATION  MEASUREMENT. • I 
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CO  TO  1020 
1010  tlRlTEO.lOlSI 

1015  FCR»*ATl//3X.*TCS/SQ.W.l..  OF  Rl- 1  POLARIZATION  ME  ASUREMENT«  ■  I 
lOZO  RETURN 
ENC 
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SUBROUTINE  MEASURICIM.  THETA  .PHI  .OilTNOO.  KCONN.NBOUNO.NNOOES. 
SHECCES.NFACES.KUHKKS. ITRAK  I 

C  THIS  SUBROUTINE  COMPUTE  THE  MEASUREMENT  CURRENT  VECTOR* 

IMPLICIT  COMPLEX  lO 
REAL  COS .CABS 

COMPLEX  CIMINUNKNS) •HMTaHMP 
COMPLEX  HX*HV*H2*HOaTT.CS(3l 
DIMENSION  OATNCC(NNOOES.3|*TMATO*2I 
INTEGER  NCCNN(NEOGES*3l •NBCUNOC 150.«)*ITRAKCNE0GES) 
COMMON/KKK/AK .P I 
CCMMCN/POLARM/HMT.HMP 
COMMQN/OIELEC/OIA.OIB 
ClaCMPLXI  1. 0*0*0 
C 

DO  2  IslvNUNKNS 
CIMd  |sCMPLX(0*0*C*0) 

2  CONTINUE 

CTIsCMPLXf cost THETA). 0*01 
CT2sCMPLX(SIN( THETA ) *0*0) 

CPHllsCMPLXCCCSIPHI 1*0*0) 

CPH12sCMPLX(SlN|PHl )«0*0) 

C 

HX*HMTMCT1 ACPHI 1-HMPACPHI2 
HVaHMTMCT 1 *CPH I2HHMPACPHI 1 
H2«-HMT*CT2 
NSEsNEOGES-NUNKNS 

c 

C  NOH  CALCULATE  THE  PARAMETERS  OF  THE  FIELD  TRIANGLE* 

DO  499  IJ«1*NFACE£ 

C  OBTAIN  THE  EDGES  OF  THE  FIELD  TRIANGLE* 

j2»necunottj»ei 

J3«NeOUNO{ I J*3) 

JAsNBOUNDI 1J*4I 

C  OBTAIN  THE  VERTICES  OF  THE  FIELD  TRIANGLE* 

IFINCONNI J2*2)*Ea*NC0NN(J3*2))  GO  TC  250 

IF(NCONN(J2*2).Ea*NCONNI J3*3> )  GO  TC  250 

NJlsNC0NN(J2.3) 

GO  TO  29S 

250  NJI^NCONNI J2*2) 

255  IFINCONNC J3*2)*Ea*NCaNN<J4*2)l  GO  TO  250 

lF(NCONN(J3.2).Ea.NCONN( J4*3) )  CO  TC  2S« 

NJ2sNC0NN<  J3.3 I 

GO  TC  258 

256  NJ2sNCCNNlJ3*2) 

256  IFINCQNNI J4*2)*Ea*NCaNN( J2*2))  GO  TO  259 

IF(NCCNN(J4*2).EO*NCONN(J2.3I)  GO  TO  259 

NJSsNCONNI J4.3) 

GO  TC  260 

259  NJ3=NCCNN<J4*2) 

260  CONTINUE 

C  OBTAIN  THE  CENTRCIO  CF  THE  FIELD  TRIANGLE* 

Xs<DATNOOINJt*2)«OATNOO(NJ2*2l«OATNOO<NJ3*2) )/3*0 
V«(DATNOOCNJl*3)HCATNaC(NJ2*3l«OATNOO(NJ3*3.)>/3«0 
C  CALCULATE  COMPCNENTS  CF  THE  TESTING  VECTOR 
C  CORRESPONDING  TC  EACH  SIDE* 

TMATII. I)«<OATNCC<NJ2*2)4CATNOOfNj3*2>)/2*0-X 
TMATII •2)s<CATNOO{NJ2*3 )40ATNC0INJ3*3) )/2*0-V 
TMAT(2« I )s(OATNCOCNJ3*2)»CATNCO|NJl*2) )/2*0-X 
TMAT(2*2)«IOATNOGlNJ3.3)HCATNCOfNJl*3))/2*0-V 
TMAT(3* t l*l0ArNC0IN31*2)«0ATNC0INJ2*2))/2*0>X 


lMAT(3«2>s(CAThCC(NJl .3 ) ♦O^TKCC IN J2 .3 n/2*0-Y 
C  OBTAIN  THE  COORDINATES  OF  THE  FIELD  TRIANGLE. 
XJt=OATKOO(NJ1.2) 

VJlsOATNOO<NJ1.3) 

XJ230ATN00(NJ2.2I 

VJ2«0ATNC0INJ2.3) 

X JSsOATNOOf NJ3.2) 

YJ3-OATNOO(NJ3.3 ) 

C  OBTAIN  THE  LENGTH  CF  THE  EACH  SIDE* 

CS( I  )«CMPLX<S0RT(<XJ2>XJ3I*«2«IVJ2-YJ3>**2).0.0> 
CS(2>sCRPLX{S0RTCCXJ3-XJl >«*2«l VJ3-YJ1>«*2).0.0) 
CSf3)sCXPLX(SCRTf f X Jt-XJ2>**2«lYJl-YJ2>*«2).0*0) 

C  COMPUTE  THE  TESTING  TECTQR  ASSOCIATED  hITH  EACH  EDGE. 

C  AGAIN  IF  ANY  OF  THE  ECCES  OF  FIELD  TRIANGLE  IS  A  BOUNDARY 
C  EDGE  THEN  THE  POINTER  JUMPS  OUT  OF  THE  LCCP. 

DO  4S0  IRsl.3 
IFfIR.EC.il  J1SJ4 
IF<IR.EQ.2)  J1«J2 
IFIIR.EG.S}  J13J3 
L  1*0 

IFCNSE.EO.Ol  GC  TC  405 
DO  390  Jsl.NSE 
IFIJt.EQ.lTRAKlJl I  GO  TO  450 
390  CONTINUE 

DO  399  Ksl.NSE 

IFIJI.GT.ITRAKIKI I  GO  TO  39? 

GC  TC  405 
397  LlaLl«l 
399  CONTINUE 

405  CTI«CMPLX(TMAT(IR«t I.O.OI 
CT2«CMPLX<TMAT(IR.2}.0.0I 
C 

ARGMNT8X*SIN<THETAI*COS<PHIl«V*SINf THETAItSINfPHI) 

CARG«CMPLX<O.0«-AK«ARGMNTI 

H00TT=HX4CT 14HY4CT2 

CHTEMPaHOCTT«CEXPfCARGI 

IFIlR.EO.l)  GO  TC  420 

IFIIR.EG.BI  GO  TO  430 

CFLAGsCI 

IFINCONNf JI.2).E0.NJ2.AN0«NCaNN<Jt.3I.EC.NJll  CFLAGs-CI 
GO  TO  440 
420  CFLAGsCl 

IFCNCONNI J1.2I.Ea.NJ3.AN0.NCCNNf J1.3).EC.NJ2>  CFLAC*-CI 
GO  TO  440 
430  CFLA6«C1 

IFfNC0NNIJ1.2I.Ea.NJl.AN0.NCCNNCJ1.3I.EC.NJ3l  CFLAG«-CI 
440  CONTINUE 

CIMf Jl>LI|sClM( Jl-Ll >4CFLAG4CS€IR>*CHTEMP 
450  CONTINUE 
460  CONTINUE 
499  CONTINUE 

999  CCNTINUE 

DO  1000  Ist.NUNKNS 

1000  CIMf 1 )«CIMf 1>*CXPLX(2.0.0.0) 

RETURN 

END 


SUeRQUTl^E  SCAINT<Xt.Yt«X2.V2«X3*V3«X«V.CPHl«AREA| 

C  THIS  SUBROUTINE. MI THE  HELP  CF  SUBRCUTIAE  lATCPL. 

C  EVALUATES  THE  SCALAR  POTENTIAL  INTEGRAL  OVER  A 
C  TRIANGULAR  REGION.  FOR  OET A ILS. PLEASE  REFER  TO  THE  NOTE. 
IMPLICIT  CCMPLEX  CO 
REAL  CABS. COS 
CCMMCN/KKK/AK. P I 
CCMMCN/VEC/XSICTl.ETACT) 

XSI<1)=1. 0/3.0 
XSIC2)>0.0SS7IS67 
XSI <3 1*0.4701 4200 
XSIC4)*XSI(3> 

XSI<S}s0.707426SS 

XSICei*0.10l26CSI 

xsif7i*xsi(ei 

ETAII IsXSICl  I 
£TAC2)3XSIC31 
ETAI3ISXS1I2I 
ETAI4)*XSIC4> 

ETAC  S}*XSI  CO 
ETAC6>=:XSI(S} 

ETAC7I*XSIC7) 

CFkCMPLXCO.O.O.OI 
00  120  1^1.7 

RI*C<X-X1|-CX2-X1 I4XSIC I 1-CX3-X1I4ETAC 1)14*2 
R2>CCV-V1)>( V2-V1 ITXSICI )-CV3-Vt )*ETACI ))**2 
R*SaRT<Rl4R2) 

CR*CMPLXC0.a.>1.04AK*R) 

IFCCABSCCR).LE.1.0E<-00  GO  TO  102 
CFl*CCEXPCCR)-CMPLXCl.O.O.aiJ/CMPLXCR.O.O) 

60  TO  103 

102  CFI*CMPLX(0.0.>AKI 

103  IFCl.EO.I)  GC  TC  lOS 

IFC l•Ea•^.aR«l•Ea•3.OR•l•E0•4>  GO  TO  110 
CF«CF4CF14CMPLX(. 1250392. 0.0) 

GO  TC  120 

105  CF3CF4CF14CMPLXC0.22S.0.0) 

CO  TC  120 

110  CF«CF»CF1*CMPLXC. 1323942.0.0) 

120  CONTINUE 

CALL  INTGRLCX1.V).X2.V2.X3.V3.X.V.PCT.AREA) 
CPHI«CF*CMPLXCAREA.O.O)TCMPLXCPOT.O.O) 

150  CONTINUE 
RETURN 
ENO 
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SUBROUTIKE  vEcihT<xi .VI .xa.vs.xa.vs* 

«X« Y.CAXSl.CAETA.AREA) 

C  THIS  SUeROUTlHE.W  ITH  THE  HELP  Of  SUBROUTINE  LIMNT. 

C  EVALUATES  THE  VECTCR  POTENTIAL  INTECRALS  EVER  A 
C  TRIAN6ULAR  REGION.FCR  OETA ILS .PLEASE  REFER  TO  THE  NOTE. 
IMPLICIT  CCMPLEX  (C) 

REAL  CABS.CCS 
COMMOHyKKK/AK  tP I 
CCMMON/VEC/XSI<7).ETA(7l 
CFsCMPLXCO.0.0.0 ) 

CG^sCMPLXlO.O.O.O) 

DC  120  lsl.7 

Rls( (X-Xt)-f X2-Xt}*XSI< I)-<X3-Xt l«ETA<I )l««2' 

R2s< I V-Vl J-CV2-V1IAXSI (I )-IV3-Vl}*ETA( 1 1 »**2 
RsS0RT(RI«R2l 
CRsCMPLXIO.O.-l .0*AK*R1 
IFICABS(CR}.LE.1.0E-0«}  GC  TC  102 
CAsICEXPICR  }-CMPLX<  t.O.O.OD/CMPLXiR.O.O 
CFt«CMPLXfXSI< I ).0.01*CA 
C61»CMPLX<ETA<I )«0.0)«CA 
GO  TO  103 

102  CFlsCMPLXIO.O.-AKAXSllf >) 

CGlsCMPLXIO.O.-AKAETAf III 

103  IFll.EQ.ll  GC  TC  105 
IFII.E0.2.CR.I«Ea.3.0R«I.EC.4t  GO  TO  110 
CFSCF4CF 1  AC MPL > I • 1 259392 « 0. 0 1 
CGsCGACGlACMPLXI .1259392. 0«0I 

GC  TO  120 

105  CFsCFA-CF1*CMPLX<0.22S.0.0} 

CC3CG«CG1*CMPLX(0.225.0.0I 
GO  TO  120 

110  CF=:CFI-CF1«CMPLX<.1323942.0«0I 

C6sCG*CGl*CMPLX I .1323942.0.01 
120  CCNTINUE 

CALL  LININTIXl •V1.X2.V2.X3.V3.X.V.PCTXSI.PCTCTA.AREAI 

CAXSI«CF*CMPLX|AREA.O.O)ACMPLX|POTXSI.O.OI 

CAETAsCG*CMPLX<AREA.O.OI«CMPLXCPOTCTA.O.O> 

150  CONTINUE 
RETURN 
ENO 
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SU8RCUTINE  L  1  MKTCXt  •  VI  .XE .VSaXa.Va cXaV • 

IPOTXSl.POTEfA.AREAl 

C  THIS  SUBROUTINEattlTI-  THE  HELP  OF  SUBROUTINE  INTGRLsEVALUATES 
C  XSI/R  ANO  ETA/R  INTEGRALS  CVER  A  TRIANGULAR  REGION.  THE 
C  QUANTITIES  DEFINED  HERE  ARE  SAME  AS  THOSE  USED  IN  THE 
C  REFERENCE  SITED  IN  THE  NOTE. 

COMMON/POTEN/PCTI 

Asixz-xiiAAE+i va-vii**a 

Ba|X3-XI  )*Aa*<Y3-Vl lA^a 

Cs»-2.0*t  I  X>X1)«(X2-X1  )»(V-Vl)«IV2»Vt  )  ) 

Ds>2.0*If X-X1)*(X3-XI|H<V-V1>*CV3-V1I» 
Es2.0«(<X2-Xt}*(X3>Xl}A<V2>VII*fV3-Vtn 
Fs(X>Xl>**2«IV-Vl lAFa 

Ats(2.0«B-C«0-El«SCRT(B«OFF>«C2.0*AHC-0>E>«SORT(AHCHF> 

A2s4.0*IA«e-E} 

A3SA1/A2 

AAsA.O*|A»C>«(BHOHF}H4.0<*F*tB-C>E>-IC«O«E>*«2 

ASse.04SQRT((A«a-CI«43l 

A6aA4/AS 

IFCABSIABI.LE. l.OE-041  GO  TO  S 
ALls2.04SQRT( AAB-EIASQRTIBHOHF) 

AL2a2.0*S0RTf A46-E |4S0RT(A4C«F| 

AL3s2.a48-C40-E 
AL  4^2 .044  40-0-2 

A J1SA34A64ALOG I ABS ( < AL 14AL3 I/I AL2-AL4 ) 1 1 
AjasA34A««AL0G<ABS( (AL24AL4>/CALI-AL3II I 
CO  TO  6 

5  AJ1SA3 
AJ33A3 

6  B1SS0RTCA4C4F J 
B2«<<2.04A4Cl*Bl-C«SaRTIF||/I4.04A I 
ANUMaABSC2.0*SCRT(A)*B142.04A4C> 

0ENaABS<2.C4SaRTCA«Fl4CI 
IF<ANUM.LE.t.0E-04l  GO  TO  ID 
IF<0CN.LE.I.0E-04I  GO  TO  10 

B3«ABS<<2.0«SGRT< A)«8142.04A4CI/f2.0*SORTI A4FI4CII 
A63sALaG<83l 

AJ43B24I4.04A4F-C442I4AB3/C6.04SQRTCA443I) 

60  TO  It 

10  AJ43B2 

11  B4sSQRTI8404F1 

B8«l 12. 04B40I4B4-04SORT IF  11/14.046  I 
ANUMaA8SI2.04S0RTlB)4B442.04B4O> 

OEN«ABS<2.04SaRT<E4F>4CI 
lFCANUM.LE.I.OE-041  GO  TO  IS 
IFIOEN.LE. l.OE-041  GO  TO  IS 

eO^ABSI <2. 04S0RT<6 >48442.046401/1 a«04SORTIB4F}4DI I 
ABOsALOGieO) 

AJ2«BS4I4.04B4F-0442)4A66/<B.04SaRTI6443l> 

GO  TC  to 
IS  A32«8S 
10  CONTINUE 

POTsPOT I/C  a.04AREA I 
AR1«2.0484CAJ1-AJ2>-E4(AJ3-AJ4I 

AR23f2.04ARt-I2.04B4C-E40>4POT>/C4.04A4B-E442| 

POTXS1*2.04AREA4AR2 

AR3«4.04A4(AJ3-AJ4)-a.04E4<AJt-Aja)-t2.e4A40-E4C>4PdT 

POTETA«<2.04AREA4AR3l/<4.04A4e-E442| 

RETURN 

END 


SUBROUTINE  INTGRLIXAl  •¥«  1  •X«2«V/t2«)i«3*Y43« 
SXaV.POT.AREAl 

C  THIS  SUBRGUTINE.WITH  THE  HELP  CP  SUBROUTINE  CA.EVALUATES 
C  THE  »/R  INTEGRAL  OVER  A  TRIANGULAR  REGICN. 
COMMON/PCTEN/PCT 1 
FI»2.0AARSINf 1.01 
XlsXAl 
Tl*VAl 
X2aXA2 
V2*VA2 
X3aXA3 
V3ayA3 

5  AR3s(X2-Xll*f V3-yt I-<V2-V1I*(X3»X1> 

AREA»ABS(AP3 1/2.0 
UN2«AR3/f2.0*AREA) 

xo»x  . . 

TO=aV 

RMsSORTMXl-XC  >**2*<V1-Y0  1A*2> 

IF(RM.6T.1.0E-06}  GO  TO  12 

X0UPPVSX2 

V0UPMYSY2 

X2«X1 

Y2*YI 

X13X3 

Yl=sY3 

X3SX0UMRV 

V3SV0UMMV 

GO  TO  S 

12  URX=IXl-XOI/RR 

ORY=< Yl-VOl/RR 
UTYs-UNZAURV 
UTVsUNZ*URX 

XTI«<X1-X0>*URXA<Y1-V0)AURV 

YT1»<X1-X0 lAUTX* <Y1-V0 lAUTV 

XT2s(X2-X01«URX«(Y2>Y01*URY 

YT2=(X2-X0>AUTXA<Y2-V01*UTY 

XT3«<X3-X01*URX*(Y3-V01*URV 

YT3*< X3-X0»AUTXA< Y3-V0I*LTY 

OETRR«2.0«AREA 

XSIa  <  XT3AYT 1-XT lAYTS l/DETRR 

ETA*I XTIAYT2-XT2AYT1I/0ETRN 

ZETAsl.O-XS  1-ETA 

SI0E1«S0RTC<XT2-XTII*A2H<YT2-YT1 )**2) 

SIOEZbSORTI <XT3-XT2I*A2«(VT3-VT21**21 

S I0E3«SQRT( f  XT 1-XT3 ••♦ZA lYT I-VT3  >A*2  > 

Te»»P®|XT2-XTl)*<XT3-XTl|P<VT2-YTl>*<YT3-VTll 

ANGLE IsARCOSI TENP/I SIDE 1«S10E3> 1 

TENP*f XT3-XT2)AIXTI-XT21A€VT3-YT2I*IYTI-VT2) 

ANGLEZsARCOSI TEPP/( S10E2ASI0E III 

ANGLE3«PI-ANGLEI-ANGLE2 

ERl«1.0E-06 

FLAG«0.0 

A00«ABS<XSIIAA8S<ETA)yABSCZETAI 
IFIAOO.GT.Il.OAERIl)  GO  TO  SO 

IFIXSI.GE.d.O-ERlI.ANO.XSl.LE.f  l.OAERI  I)  GO  TO  IS 
IFtETA.GE.d.O-ERII.ANO.ETA.LE.Cl.O^ERllI  GO  TO  20 
IF(ZETA.CE.d.O-ERll.ANO.ZETA.LE.(I.O«ERll}  GO  TO  2S 
IF<XSI.GE.-ER1.AN0.XSI.LE.ERII  6C  TC  30 
IFIETA.GE.-ERI.ANO.ETA.LE.ERII  GO  TO  35 
IF<ZETA.GE.-ERl.ANO.ZETA*LE.ERtl  GO  TC  40 


FLAGsl.O 
€0  TC  SO 

IS  CALL  CAIXT3.VT3.XTI.VTI.VAL1> 

VAL=VAL1 
CO  TO  too 

20  CALL  CA(XTl.VTl.XT2.VT2.VALli 

VAL=VAL1 
CO  TO  IOC 

2S  CALL  CA<XT2.VT2.xra«vr3.VALt> 

VALsVALl 
CO  TO  too 

30  CALL  CA< XTt.V11.XT2.VT2.VALI) 

CALL  CA(XT2.VT2.XT3.VT3.VAL2> 
VAL3VAL1AVAL2 
GO  TO  too 

35  CALL  CAIXT2.VT2.XT3.VT3.VAL1) 

CALL  CA<XT3.VT3.XTt.VTt.VAL2) 
VAL»VALt  4VAL2 
GO  TO  too 

40  CALL  CA(XT3.VT3.XT1.VT1.VAL1) 

CALL  CAIXT1.VT1.XT2.VT2.VAL2I 
VAL«VAL1 4VAL2 
GO  TO  too 

50  CALL  CA(XTl.VTt«XT2«VT2.VALl) 

CALL  CA(XT2.VT2.XT3,VT3.VAL2I 
CALL  CAf XT3.VT3.XT1 .VTl.VALal 
VAL3VAL14VAL24VAL3 

too  CONTINUE 
POTsVAL 
POT l»POT 
PETUPN 
END 


SUBROUTINE  CACXl*Yl*Ha»t2«VAL» 
COMNON/ERRCR/ERR 1 
RAaSORTI Xl**2*Tt**2) 
R8*S0BT<X2**2*V2**2> 

iH.=S0RT<<X2-XI  l♦♦2♦<r^-V»I♦•2> 

OOTs(  *-XI>*{  X2-XlI*<->rII*<V2-Vl  >>XXL 
X0sXt«00T*(X2-Xll/AL 
V0=Vl*D0T*«T2-Tl >/AL 
RN0TaSaRT<X0*A2«Y0**2) 

ERRIaRNCT/AL 

2eROaI-0E-06 

IE(ERRl«L£«2ERO)  GO  TO  tC 
PHII«ATAA2<Yt.Xl ) 

PH12aATAN2(V2«X2) 

RHlNOTaATAN2< VOaXOl 
FtaEXPRNCRNOT.PHiNOT.PNill 
F2aEXPRN<RNaT«PHlNaT«PHI2 I 
VAL=F2-F1 
cc  TC  n 

10  VALsO«0 

1 1  RETURN 
ENO 
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FUNCTION  EXPRMFNCT.PHlNOTvPMII 
COMMON/ERROR/ERP 1 
TEPFIaSCRT (RN0T*«2> 

TEPP2«RN0T*SlMPHIN0T-PFi  I 
IFCERRI.LE.  t•0E->06l  GO  TO  10 
/ILPHA|3PH1N0T<-P»i1 
ALPHA2sARSIN< 1.0} 

ERR2S AES ( ALPHA 1**2-ALPFA2««21 
IF(ERR2.LE«I.0C-061  60  TC  10 
TEMP3«ALC6« ( TEMP  1 *TEMP21/(TEPP1>TEMP2 ) ) 
CC  TC  11 

10  TEMP3S0.0 

11  EXPRNs-<RhCT*TEPP3)/2.0 

RETURN 
END 


V 
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SUBROUTINE  CSM INV C A«NOIM» N»CETERM«CONO« lERRl 

COMPLEX  Af KOI M.NOIM} .PIVOT  1250 l•AMAX•T•SMAP•OETERM•U•CMPLX•CONJG 

1NTEGERA4  IP I VCT ( 250 1 . I KOE X 1250 .2 1 

REAL  TEMP.ALPKACBSOl.CABS 

COMPLEX  CTEMP.CALPHA(250 > 

lERRsO 

IFIN0IM.LE.2501  GC  TO  5  ^ 

I ERRS 1 

MRITEO.AI  KOIM 

4  FORMAT! •OCSMINV  ERROR.  ATTEMPT  TO  INVERT  A  MATRIX  ■ 14. 

!•  CK  A  SIDE.*/*  MFEN  250  X  250  IS  ThE  MAXIMUM  ALLGMEO.*) 

RETURN 

5  CONTINUE 

OETERM  s  CMPLXIl. 0.0.01 

SUMAXAsQ. 

CO  20  JSI.N 
ALPHA<J)=0.0 
CALPHAI J)sf O.O.C.O ) 

SUMROMsO. 

OC  10  lal.K 

CALPHA!  J)sCALPHA(  JiAAf  J.D*  CCNJGIACJ.  ID 
ALPHAIJisREALCCALPHAI J1 ) 

10  SUMRChsSUMPOII  ♦  CABS(A(J.1}| 

ALPHAIJis  SORTIALPHAI J)  } 

IFISUMROH.GT.SUMAXAI  SUMAXAsSUMRCM 
20  IPIVCT(J)=0 

00  600  1=1. N 
AMAX=CMPLX(O.O.C.O) 


OO 

105  J=1.N 

IF 

IIPIVOT! J)-ll 

60. 

105. 

60 

00 

100  K=1«K 

IF 

!1PIV0T!K)-11 

80. 

too. 

740 

80  CTEMP=AMAX«  CON  JG  I  AMAX  )>A<  J.K  }«  CONJGCACJ.KM 
TEMPsREALICTEMPI 
lFfTEMP18S.es. ICO 
85  IRC)l=J 

ICOLUMsK 
AMAX=A(J.K  I 
100  COKTIKUE 
lOS  CONTINUE 

IPlVOTf ICOLUMIslPIVOTflCOLLMltl 
IF  flROM-ICOLUMI  140.  260.  140 
140  OETERM=-OETERM 
OO  200  L=1.N 
SMAP=A( IRCM.L} 

AIlROR.LlsAflCCLUM.Ll 
200  A<  ICOLUM.L  )  =  SllAP 
SMAPsALPHAf  IROV I 
ALPHA! IRO«)=ALPNAC ICOLUMl 
CALPHA! 1C0LUM1=S«AP 

ALPHA! ICOLUM 1=RE AL  fCALPHA ! ICOLUMl 1 
260  INDEX!  1.1  MIROM 

INDEX! 1.2I«1C0LLM 
PIVOT! I }=A! ICOLUM. ICOLUMl 
U  «  PIVOTIl) 

ALPHAlsALPHAIlCCLLMI 
CALL  DTRMNTfOETERM.U.ALPHAII 
CTEMP=P1V0T!1 M  CONJGIPIVOT! 1 1 1 
TEMPsREALICTEMP) 

IF !TEMP 1330.720.330 


r 
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330 

35C 

380 

400 


4S0 

550 

600 

620 


630 


705 

710 


900 

910 


720 

730 

740 


A( ICOLUM.ICCLUV)  3  CMPLXl l.OaO.Ol 
DC  350  L=1«K 
U  s  PIVCTll) 

AdCOLUM.L)  =  Af  ICCLUP.Ll/t 
CO  550  L1=1»N 

IFCLl-ICOLUM)  400*  550.  400 
T«A(1.1»1C0L(;P) 

AILI*  ICOLUP)=»  CPPLX(0»0.0«0) 

CC  450  L=1*K 
U  «  AdCOLLM.L) 

A<L1*L}  *  ACLI.1.1-U4T 

CCKTIKUE 

COKTINUE 

00  710  Isl.N 

L*K41-I 

IF  dNOEXIL,!  )-IK0eX<L.2» )  630«  710.  «3C 
JR0«=:IN0EX<L»1  ) 

JCCLUP= IN0EX(L«21 
00  70S  K=1.N 
ShAP^AIK. JPC4) 

A<X«  JPOttlsAIK.  JCCLUM) 

AIK. JCCLUM1«SHAP 
COKTINUE 
CCKTIKUE 
SUKAXI^O. 

00  910  Isl.K 
SUMPOV=0. 

CC  900  J«1.N 

SUMPOksSUMPOli  *  CAeSIAd.Jl) 

IFISUMftOW.CT.SUPAXl  1  SUMAXI*SUMRC» 

COKTINUE 

CONO  -  1./<SUKA>A*5UMAX11 
RETURN 

FORNATl’o'dOl  ••♦•♦♦♦♦•1/*0RATR1X  IS  S I  K6UI.AR* /‘O*  •  10  I  > 

RETURN 

ENC 
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SUBROUTINE  D  f  RRN  T  <  OETEPN  .U  .il  » 

real  cabs 

CCRPLEX  DETERM.U.CNPLX 
CORWON/SCAFAC/ISCALE 
DATA  lSCALE/0/ 

IFICABSIOETERM)  -GT,  t.E~10)  GO  TC  100 

OETERM=OETERM« 1 .E  t 0 

ISCALE-ISCALE^l 

100  OETERMsOETERR«U/CRPLX<A.O.OI 
return 

ENO 


AD-AUb  593 


UNCLASSIFIED 
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